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ABSTRACT 

•J*;  I. 

Molecular  dynamics  are  computed  for  model  atom  transfers  A  +  BC  AB  +  C 
in  rare  gas  solvents  at  liquid  densities.  We  find  that  the  reaction  dynamics  can  be 
understood  in  terms  of  a  simple  picture  which  consists  of  three  stages:  1)  activation  of 
reactants,  2)  barrier  crossing,  and  3)  deactivation  of  products.  The  effects  seen  in 
stages  1)  and  3)  can  be  largely  interpreted  in  terms  of  existing  models  of  energy  and 
phase  decay  in  solution,  while  the  effects  seen  in  stage  2)  can  be  largely  interpreted  in 
terms  of  gas  phase  A  +  BC  barrier  crossing  dynamics.  We  find  that  Transition  State 
Theory  is  in  perfect  agreement  with  the  simulations  for  the  20  and  10  kcal/mol  barrier 
reactions  and  is  a  very  good  description  for  a  S  kcal/mol  reaction  barrier.  At  low  bar¬ 
rier  curvature  dynamical  effects  due  to  the  solvent  are  shown  to  induce  some  recross¬ 
ings  of  the  transition  state  barrier,  thus  causing  rate  constants  calculated  by  simple 
transition  state  theory  to  be  slightly  too  high.  A  modification  of  transition  state 
theory,  which  considers  the  effect  of  the  time  dependent  friction  of  the  solvent  on  the 
dynamics  at  the  transition  state,  is  shown  to  predict  corrections  to  the  rate  constants  in 
very  good  agreement  with  the  results  from  the  simulations,  ip  4  i  >  i  , «  r  \  -  "D 
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I.  Introduction 

The  simplest  chemical  reaction  involving  the  making  and  breaking  of  different  chemical  bonds  is 
the  triatomic  atom  transfer  reaction  A  +  BC  -»  AB  +  C.  Its  experimental  and  theoretical  study1  has 
been  central  to  the  detailed  understanding  of  gas  phase2'7  reactions,  contributing  for  example  to  the 
development  of  Transition  State  Theory  (TST)1>8  -  which  gives  in  many  cases  a  clear  and  accurate  pic¬ 
ture  of  gas  phase  reaction  kinetics  -  and  to  the  development  of  Molecular  Dynamics  (MD)  to  under¬ 
stand  the  detailed  microscopic  motions  involved  in  reactions.  While  the  theoretical  understanding  of 
the  microscopic  processes  of  gas  phase  reactions  is  now  well  advanced,  the  corresponding  understand¬ 
ing  of  the  detailed  aspects  of  reactions  in  liquid  solution  is  still  sketchy,  in  spite  of  what  is  arguably  the 
greater  relative  importance  of  solution  reactions  in  chemistry.  (For  reviews,  see  Refs.  8,10-lS.) 

The  present  work  is,  to  our  knowledge,  the  first  full  scale  molecular  dynamics  study9  of  the  A  + 
BC  reaction  in  solution  in  which  reactive  trajectories  are  computed.  We  analyze  the  resulting  classical 
trajectories  to  try  to  answer  some  important  questions  about  the  nature  of  the  A  +  BC  reaction  in  solu¬ 
tion,  including: 

1)  What  is  the  character  of  ABC  trajectories  that  lead  to  reaction,  and  what  is  die  solvent 
influence  thereon?  We  find  that  these  trajectories,  except  at  the  barrier  top,  are  quite  different  in  the 
solution  and  gas  phases  and,  away  from  the  barrier,  can  be  characterized  in  terms  of  energy  transfer  and 
phase  decay  characteristics  of  various  (eg.  translational,  rotational,  vibrational)  degrees  of  freedom  of 
the  nascent  reactant  and  product  species.  In  addition,  the  role  of  the  solvent  momenta  in  the  reaction 
process  is  examined  through  solvent  momentum  perturbation  calculations. 

2)  When  is  TST  valid*"15  in  solution?  The  most  apparently  vulnerable  feature  of  simple  TST  is 
the  assumption  that  trajectories  originating  from  reactants  and  achieving  activation  energy  will  always 
go  on  directly  to  form  products  with  no  recrossings  of  the  transition  state.4  Clearly,  this  assumption  is 
in  jeopardy  if  the  situation  arises  in  which  a  strongly-interacting  solvent  impedes  successful  barrier 
crossing  through  "collisions"  with  the  reacting  complex  near  the  transition  state,  thereby  inducing 
recrossings  and  reducing  the  rate  constant  below  the  TST  value.  One  might  term  this  a  dynamic  solvent 
"cage"  effect  in  the  transition  state  neighborhood.  We  examine  this  question  in  some  detail,  and 
analyze  our  results  in  terms  of  the  Grote-Hynes  perspective  of  a  time  dependent  friction16"18  that  the 
solvent  exerts  on  the  reaction  coordinate  in  the  neighborhood  of  the  transition  state,  and  its  effect  on  the 
rate  constant. 

The  appealing  feature  of  TST  is  the  simple  but  intuitive  picture  it  gives  for  understanding  tire 
nature  of  reactive  barrier  passage:  an  activated  complex,  in  quasi-equilibrium  with  reactants,  proceeds 
directly  to  farm  products.  A  critical  assumption  of  TST  -  the  existence  of  an  equilibrium  between 
reactants  and  activated  complexes  -  can  be  realized  in  a  solvent  where  continual  collisions  occur  with 
the  reacting  system,  although  these  same  interactions  might  sometimes  cause  the  no-recrossing 


ity  Codes 

and/or 

social 


-2- 


condition  to  break  down.  If  a  trajectory  recrosses  the  transition  state  surface  then  TST  will  overestimate 
the  rate  constant  because  it  assumes  that  every  crossing  of  the  barrier  toward  products  contributes  to  the 
overall  rate.  This  recrossing  effect  is  accounted  for  by  the  so  called  transmission  coefficient  k  correc¬ 
tion  to  Transition  State  Theory.  As  is  well  known,  the  actual  rate  constant  k  is  thus  written  in  terms  of 
the  simple  TST  rate  constant  A7*7  as 

k  =  *nrK,  (1.1) 

where  K  adjusts  the  TST  result  to  correct  for  the  solvent  induced  recrossings.  Some  basic  ingredients 
which  determine  k  are  considered  in  a  highly  idealized  treatment  of  reactions  in  solution  by  Kramers,19 
in  which  the  reacting  system  is  modeled  as  an  effective  particle  of  mass  it  moving,  in  the  transition 
state  region,  on  an  inverted  parabolic  potential  of  imaginary  frequency  whose  real  magnitude  is  to*  and 
subject  to  solvent  friction  £.  The  latter  is  defined  in  terms  of  the  time  correlation  function  of  the  ran¬ 
dom  force  exerted  by  the  solvent  on  the  reaction  coordinate, 

m 

;  =  JdfC(r),  (1.2) 

in  which  £(r)  is  given  by  the  fluctuation-dissipation  theorem  as 

ft*)  =  Ff(‘) >  .  (1.3) 

where  F  is  the  solvent  force  on  the  effective  particle  of  mass  p.  and  (  )  indicates  a  solvent  phase-space 
average,  kB  is  Boltzmann’s  constant  and  T  is  the  temperature.  The  Kramers  theory  emphasizes  the  key 
parameter  £  /  to4  in  the  determination  of  k.  If  the  frictional  forces  are  weak  compared  to  the  intrinsic 
reaction  forces  acting  on  the  particle  in  the  barrier  region,  then  £  /  cot  «  1,  k  -»  1  and  k  -*  A757.  If, 
on  the  other  hand,  the  barrier  curvature  is  very  small  and  die  retarding  influence  of  the  solvent  is  large, 
then  £  /  to*  »  1  and  k  -*  to*  /  £.  In  this  limit  the  solvent  is  very  effective  in  inducing  recrossings  and 
the  passage  from  reactants  to  products  over  the  barrier  is  essentially  diffusion-controlled. 

In  die  latter  case,  and  in  regimes  approaching  it,  the  reaction  rate  is  affected  by  the  long  time 
macroscopic  friction  in  the  Kramers  description.  Than  are,  however,  important  cases  for  which  this 
simple  picture  due  to  Kramers  breaks  down.  It  has  been  argued  by  Grote  and  Hynes16'18  that  TST  can 
in  fact  be  a  good  approximation,  even  for  strong  friction,  if  the  barrier  curvature  is  very  sharp.  For  this 
case  the  reacting  system  spends  so  little  time  on  the  barrier  (on  the  order  of  time  o£l)  that  the  solvent 
is  unable  to  respond  to,  and  retard  motion  across  the  barrier  since  there  is  not  enough  time  for  any 
effective  collisions  to  occur.  The  effective  solvent  friction  at  the  transition  state  is  small  because  the 
short  time  scale  solvent  response,  given  by  the  early  part  of  the  time  dependent  friction,  £(r),  is  more 
important18  than  the  long-time  overall  response,  £.  We  will  see  that  this  is  an  especially  important  per¬ 
spective  for  sharp  barriers. 

The  solvent  can  also  play  another  important  and  quite  distinct  role;  dissipating  the  excess  energy 
of  newly  formed  products.  If  energy  deactivation  of  products  is  achieved  a  short  time  after  the  transi¬ 
tion  state  is  crossed,  then  it  is  unlikely  that  there  will  remain  sufficient  energy  to  recross  the  barrier, 
this  helps  TST  to  correcdy  describe  the  reaction.  Conversely,  the  lack  of  sufficiently  rapid  energy 
transfer  will  lead  to  barrier  recrossing  by  the  energetically  "hot"  nascent  products.  In  this  event,  the  rate 
constant  would  depend  on  die  rate  of  energy  stabilization  of  the  nascent  products  by  the  solvent 

Similar  to  the  role  that  the  solvent  plays  in  the  energy  deactivation  of  products  is  the  role  it  plays 
in  providing  energy  for  the  activation  of  reactants.  It  is  conceivable,  for  example,  that  this  activation 
might  be  slow  compared  to  the  barrier  passage  step,  so  that  activation  dynamics  explicitly  enter  the 
rate.  It  is  also  possible  that  only  a  small  amount  of  solvent  phase  space  can  be  effective  in  the  activa¬ 
tion  of  the  reactants  to  form  products. 

One  purpose  of  this  paper  is  to  investigate  the  validity  of  Transition  State  Theory  for  atom 
transfer  reactions  in  monatomic  liquids.  We  simulate  reactions  using  molecular  dynamics  and  explore  a 
range  of  the  system’s  properties  (characterized  by  £(/),  to*,  etc...),  discovering  both  where  TST  can  be 
applied  successfully  and  where  it  breaks  down.  In  addition,  we  search  for  a  simple  picture  to  understand 
the  reaction  molecular  dynamics,  and  begin  the  process  of  investigating  the  detailed  role  of  the  solvent 
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in  affecting  the  course  of  the  reaction. 

The  outline  of  the  paper  is  as  follows.  In  Sec.  II  we  present  in  detail  the  computational  techniques 
used  to  simulate  A  +  BC  trajectories  in  rare  gas  solvents.  Section  m  explores  the  reaction  X  +  X2  in 
argon  solvent  In  Sec.  IV  we  examine  the  effect  of  varying  the  nature  of  the  solvent  Concluding 
remarks  are  given  in  Sec.  V. 

n.  Computational  Techniques 

Molecular  dynamics  techniques  have  been  previously  used  in  several  condensed  phase  reaction 
simulations;  examples  include  defect  motion  in  solids,20’21  isomerization22  and  radical  recombina¬ 
tion23*24  in  solution,  transformations  in  macromolecules, 251 26  and  adsorption  on  surfaces.27  The  solution 
phase  MD  simulation  of  the  A  +  BC  reaction  requires  special  attention  to  initial  sampling  and  trajectory 
techniques  which  we  now  discuss. 

A.  Molecular  Dynamics 

Molecular  dynamics  are  used  to  compute  time  histories  of  the  positions  and  momenta  for  all  of 
the  atoms  in  the  system.  This  is  accomplished  by  numerical  integration  of  the  equations  governing  the 
classical  morions  of  particles  in  a  conservative  force  field,  i.e.  Hamilton’s  equations: 

0.1) 


in  which  pN  and  r*  are  the  conjugate  momentum  and  position  for  a  system  containing  N  unique  parti¬ 
cles,  and  H  is  the  Hamiltonian  for  the  system.  Given  a  set  of  position  coordinates,  r*  the  force  on  any 
particular  atom  is  computed  as  the  sum  of  the  forces  over  all  pairwise  additive  interactions,  and  a  3- 
body  term  is  used  to  describe  the  internal  A  +  BC  interaction  (see  below).  A  modified  Verlet 
algorithm,28'30  incorporating  a  rime  step  of  1.0  femtosecond,  is  used  to  integrate  these  equations  of 
morion,  and  truncated  octahedron  periodic  boundary  conditions31  are  used  in  order  to  approximate  an 
infinite  liquid.  There  are  100  rare  gas  solvent  atoms  used  in  the  liquid  simulations. 

These  calculations  are  carried  out  on  a  Floating  Point  Systems  AP120B  array  processor32  attached 
to  a  VAX  11/780  host  computer.  One  thousand  time  steps  require  approximately  100  seconds  of  real 
time  on  the  array  processor.  The  calculations  presented  here  required  IS  days  of  array  processor  time, 
and  would  have  required  17  months  of  VAX  1 1/780  time  if  the  array  processor  were  not  used. 

B.  Potential  Energy  Surface 

In  explicitly  constructing  the  Hamiltonian  in  Eqs.  (2.1)  and  (2.2),  a  potential  energy  surface  is 
needed  to  describe  the  interactions  between  all  of  the  atoms  in  the  system.  The  potential  energy  surface 
used  here  is  intended  to  represent  a  liquid-phase  rare-gas  solvent  in  which  the  symmetric  A  +  BC 
hypothetical  atom  transfer  reaction,  X  +  X2  — »  X2  +  X,  occurs.  The  X  atoms  are  modeled  as  uncharged 
hidogen-like  atoms  with  the  mass  of  chlorine  33.  We  stress  that  this  is  a  hypothetical  model  of  a  sym¬ 
metric  atom  transfer  (The  actual  interhalogen  transfer  reactions  appear  to  often  involve  low  barriers  and 
bound  intermediates.3 ). 

We  choose  a  surface  constructed  as  a  combination  of  hypo-surfaces  representing  smaller  parts  of 
the  problem.  First,  the  potential  energy  of  the  interaction  between  any  atom,  i,  and  any  solvent  atom,  j, 
is  represented  by  the  Lennard-Jones  6-12  potential 
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in  which  is  the  distance  between  the  atoms,  ev  is  the  depth  of  the  minimum  in  the  potential  and  ov 
locates  the  finite  mtermolecular  distance  at  which  -  0.  The  values  of  £tJ  and  ov  used  for 
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various  rare  gas33  and  X-X  interactions34'35  in  these  calculations  are  given  in  Table  1.  Table  1  gives 
parameters  only  for  the  homonuclear  pair  interactions.  For  the  heteronuclear  interactions  the  values  of 
Cy  and  Oy  are  calculated  using  the  combining  rules 

ty  =  (tg  tj,)*  ,  (2.4) 

Oy  =  (aM  +  OJ/2.  (2.5) 


Table  I.  Lennard-Jones  Parameters 

-w-fen 

Atom  -  Atom 

Oy 

eJk, 

Iatendioo 

A) 

(Jlmol) 

CK) 

He -He 

2.57 

89.79 

10.80 

Ar-  Ai 

3.40 

997.48 

120.00 

Xe-Xe 

3.85 

201120 

249.83 

X-X  (X-Ct) 

3.12 

2868.68 

345.02 

Second,  the  potential  energy  function  used  for  the  X-X-X  interaction  is  the  3-body  London- 
Eyring-Polanyi-Sato  (LEPS)5  surface 

=  Qi  +  Qi  +  Cs  -  (i?  +  A  +  A  ~  -  Jifi  ~  •Vi)'*  *  (2.6) 

in  which  r,  and  r2  are  the  variable  bond  lengths  between  the  "middle"  X  atom  and  the  two  "outside"  X 
atoms,  r3  is  the  bond  length  between  the  two  "outside"  X  atoms,  and  &  and  J,  are  linear  combinations 
of  the  "singlet"  and  "triplet"  diatomic  energies: 

=  ('£,  +  %Y2  ,  (2.7) 

and 

Hr,)  =  ('£,  -  •  (2.8) 

Here  the  "singlet"  ground  state  energy  is  represented  as  a  Morse  potential, 

%<r,)  =  'a|i  -  exp  [-‘3, A  -  'r?)]|  -  >A  , 
while  the  "triplet"  diatomic  energy  is  an  anti -Morse  potential, 

’£.<'.)  =  JD,|l  +  exppfcto  -  3r?)]|  -  3A  • 

The  LEPS  potential  is  an  analytic  form,  convenient  for  use  in  molecular  dynamics  simulations, 
which  displays  excellent  asymptotic  properties,  reducing  to  a  diatomic  bound  Morse  potential  when  one 
atom  is  infinitely  separated.  In  addition,  it  is  highly  adaptable,  containing  18  non-symmetry  adapted 
parameters  (lD„  *D„  'p„  Jp„  ‘if,  Vf,  for  i  -  1,2,3),  allowing  the  simultaneous  adjustment  of  barrier 
height  and  barrier  frequency  o>»  (i.e.  the  magnitude  of  the  imaginary  frequency  of  an  inverted  parabolic 
approximation  to  the  potential  at  the  saddle  point  for  a  linear  geometry).  We  report  all  frequencies,  <X)b, 
in  spectroscopic  units,  U.  (to*  /  2jut)  cm-1.  Table  2  shows  the  values  of  the  18  LEPS  parameters  used 
in  these  simulations,  producing  potentials  with  specific  barrier  heights  and  frequencies.  In  all  cases,  the 
Xi  diatomic  potential  in  the  limit  of  the  X  atom  far  away  is  a  Morse  potential  fitted  to  the  G2  diatomic 
potential.36  Figure  1  illustrates  potential  energy  contours  for  a  20  kcal/mol  barrier  height  and  linear 
geometry.  The  coordinate  axes  are  inclined  at  an  angle  of  60°  rather  than  90°  since,  for  the  X3  system, 
this  choice  of  axes  diagonalizes  the  kinetic  energy.1  The  equations  of  motion  of  the  reaction  system  are 


(2.9) 


(2.10) 


Table  2.  X3  LEPS  Potential  Parameters 


Barrier 

Height 

(kcal/mol) 

20 

10 

5 

Barrier 

frequency 

(cm-1) 

419 

288 

168 

(A) «  -  1,2,3 

1.9870 

1.9870 

1.9870 

M  (A)  i  «  1  A3 

1.9870 

1.9870 

%  (A"1) »  -  1,2,3 

2.0024 

2.0024 

3P,  (A"1)  i  -  1,2,3 

2.0024 

2.0024 

lDi  (kcal/mol)  i  -  1,2,3 

57.9740 

57.9740 

57.9740 

3D,  (kcal/mol)  i  -  1,2,3 

44.2290 

32.3965 

27.2053 

thus  represented  by  a  panicle  with  point  mass  |i  sliding  on  the  potential  energy  surface. 


C.  Initial  Conditions 


The  reaction  coordinate  of  X  +  X2  -»  X2  +  X  in  the  barrier  neighborhood  is  well  defined  by  the 
asymmetric-stretch  normal  mode  coordinate  of  X-X-X,  while  the  magnitude  of  the  corresponding  ima¬ 
ginary  normal  mode  frequency  is  the  barrier  frequency,  (0*.  Other  vibrational  degrees  of  freedom  for  X3 
are  described  by  the  bend  and  symmetric-stretch  normal  modes. 


If  we  attempted  to  discover  a  set  of  reactive  trajectories  by  starting  with  reactants  and  solvent  in 
equilibrium,  we  would  have  exhausted  the  capacity  of  the  fastest  available  computer  long  before  arriv¬ 
ing  at  our  answer.  Therefore  we  use  the  technique  of  Keck37-38  and  Anderson39-40  to  initialize  the  tra¬ 
jectory  from  an  equilibrium  distribution  of  all  variables  except  the  position  along  the  reaction  coordi¬ 
nate  which  is  constrained  to  the  transition  state  surface.  We  compute  the  dynamics  both  forward  and 
backward  in  time,9- 41  thus  determining  if  the  trajectory  is  reactive  or  unreactive.  The  precise  nature  of 
this  equilibrium  distribution  on  the  barrier  is  highly  dependent  on  the  details  of  the  potential  energy  sur¬ 
face  at  the  barrier  top.  Here,  the  initial  conditions  for  Xj  are  determined  by  Boltzmann-weighted  sam¬ 
pling  of  positions  and  momenta  in  the  appropriate  degrees  of  freedom  described  by  the  normal  modes. 
The  initial  conditions  are  chosen  to  optimize  the  calculation  of  the  ensemble-averaged  values  of  the 
properties  of  interest  An  ensemble  average,  (A),  is  defined  classically  as 


<A)  = 


ffdr»dp"Ae-HMi>T 
jj  dr"d^  e-"^****’7 


(2.11) 


where  r*  is  a  vector  of  the  N  Cartesian  coordinates  of  the  system,  pN  is  a  vector  of  the  momentum 
coordinates  conjugate  to  r*.  //(r*,p*)  is  the  full  system  Hamiltonian,  kt  is  Boltzmann’s  constant  and  T 
is  the  temperature.  We  select  r*  and  p*  from  a  Boltzmann  distribution  on  the  transition  state  surface. 
For  the  case  of  a  solvent  interacting  weakly  with  the  X3,  the  transition  state  can  be  conveniently  defined 
as  lying  along  die  symmetric  stretch  and  bend  coordinates  passing  through  the  lowest  lying  saddle  point 
on  the  gas  phase  LEPS  potential  surface  of  X3.  We  now  turn  to  the  details  of  the  initial  condition  sam- 
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1.  Hamiltonian  Separation 

The  reaction  system  coordinates  comprise  a  set  of  ‘fast*  degrees  of  freedom,  q,  associated  with 
die  X)  vibrational  normal  mode  coordinates  and  the  remaining  "slow"  degrees  of  freedom,  R,  charac¬ 
teristic  of  the  solvent  translational  motion  and  the  rotational  and  translational  motions  of  X3.  In  die 
spirit  of  the  familiar  Bom-Oppenheimer  approximation  for  separation  of  nuclear  and  electronic  motions, 
the  Hamiltonian  appropriate  for  our  initial  sampling  is  written  as 

//(r*.p*)  =  W(q,p,R,P)  =  tfi(q.p;R)  +  H2( R,Pfi) ,  (212) 

in  which  p  and  P  are  momenta  pertinent  to  q  and  R  respectively,  //i(q,p;R)  is  the  Hamiltonian  for  X3 
which  is  parametrically  dependent  on  die  solvent  coordinates  R,  and  W2(R,P;q)  is  the  Hamiltonian  for 
the  solvent  which  is  parametrically  dependent  on  the  average  value  of  the  X3  vibrational  coordinates,  q. 
If  necessary,  this  formalism  may  be  generalized42  to  allow  for  quantization  of  the  fast  variables  q. 

In  the  limit  of  weak  coupling  between  X3  and  the  solvent,  we  can  ignore  the  parametric  depen¬ 
dence  of  Hi  on  die  liquid  coordinates  R,  and  Hx  reduces  to  the  X3  gas  phase  Hamiltonian  //t(q,p).  The 
solvent  Hamiltonian,  H2,  is  parametrically  dependent  on  q  which  we  take  as  q^,,  the  coordinates  of  the 
saddle  point  Thus,  in  the  weak-coupling  limit,  H  becomes 

//(q,p,R,P)  =  //t(q,p)  +  H2(Rf  \qn)  .  (2.13) 

The  problems  associated  with  a  molecular  dynamics  evaluation  of  any  ensemble  average,  such  as 
Eq.  (2.11),  are  that  standard  random  sampling  techniques  give  a  poor  sampling  of  the  sensitive  variables 
q,  and  many  different  initial  conditions  must  be  sampled  for  proper  convergence.  In  addition,  random 
sampling  techniques  do  little  to  provide  physical  insight  into  the  regions  of  phase  space  which  contri¬ 
bute  most  to  the  integral  and,  therefore,  are  important  into  the  overall  "chemistry"  of  the  process.  For 
example,  the  lowest  lying  saddle  point  on  the  LEPS  potential  surface  occurs  for  linear  X-X-X.  This 
may,  at  first,  lead  one  to  believe  that  the  most  probable  configuration  for  a  reactive  trajectory  would  be 
a  col  linear  arrangement  of  X-X-X,  as  this  requires  a  minimum  of  activation  energy  to  cross  the  barrier. 
However  a  significantly  bent  configuration  is  actually  the  most  probable  arrangement  of  X-X-X;  this  is 
due  to  the  en tropic  effect  of  more  volume  in  configuration  space  at  bent  geometries  and  compensates 
for  the  increased  activation  energy  required.  In  fact,  the  probability  for  an  exactly  linear  configuration  is 
zero.  We  will  take  this  feature  into  account  in  what  follows. 

The  distinction  between  the  fast  and  slow  variables  in  the  system,  as  well  as  the  assumption  that 
there  is  weak  coupling  between  the  X3  moiety  and  the  solvent,  allows  us  to  choose  initial  conditions  for 
the  X3  fast  (t.e.  vibrational)  coordinates,  q,  independently  of  the  choice  of  initial  conditions  for  the 
remaining  slow  variables,  R,  (z'.e.  solvent  translation,  X3  translation  and  rotation).  In  addition,  classi¬ 
cally,  the  initial  conditions  for  momentum  variables  can  be  selected  independently  of  the  position  coor¬ 
dinates.  We  choose  an  integration  scheme  which  takes  advantage  of  the  parametric  dependencies  in  Eq. 
(2.13)  allowing  us  to  determine  initial  conditions  separately  for  the  X3  and  for  the  solvent,  taking  the 
parametric  degrees  of  freedom  into  account  in  an  average  way.  A  good  approximation  to  Eq.  (2.11) 
useful  for  generating  a  quadrature  scheme  thus  is 

J  dqdp  /(q.p)  e~H^yk*T  \  dRdP  A(q,p,R,P) 

1  J  y(q.p)  J  dKJT  .-‘vwvv  <M) 

We  write  the  X3  Hamiltonian  as 

«i(q.P)  =  *i(q)  +  Tj(p)  (2.15) 

and  the  solvent  Hamiltonian  as 

HifrJ/uJ  =  ®i(R;q«,)  +  T2(P) ,  (2.1 6) 

in  which  Oj(q)  is  the  gas  phase  potential  energy  (LEPS)  of  the  X3  reaction  system,  d>j(R;q^)  is  the  sol¬ 
vent  potential  energy  which  is  parametrically  dependent  on  q^,  T](p)  is  the  vibrational  kinetic  energy 
of  X3,  T2(P)  is  the  solvent  translational  kinetic  energy  plus  X3  translational  and  rotational  kinetic 
energy,  and  /(q,p)  is  the  Jacobian  which  is  yet  to  be  determined  for  transformation  from  Cartesian 


-7- 


variabies  to  the  internal  vibrational  coordinates,  q,  of  Xj.  Eq.  (2.14)  is  evaluated  using  a  product  qua¬ 
drature  in  which  the  q,  p,  R,  and  P  coordinates  are  chosen  independently:  P  by  a  Gaussian  random 
number  generator,  R  by  molecular  dynamics,  and  q  and  p  by  product  Gauss-Hermite  quadrature  in  the 
internal  coordinates  (see  below).  Note,  however,  that  the  quadrature  scheme  used  provides  a  formally 
exact  evaluation  of  Eq.  (2.11)  in  the  limit  of  an  infinite  number  of  trajectories.  The  approximation  used 
in  generating  the  quadrature  only  slows  the  convergence. 

2.  Effective  LEPS  Potential 

The  integral  in  the  ensemble  average  Eq.  (2.14)  over  the  X3  vibrational  degrees  of  freedom  can 
now  be  written  as 

J-Mp 4W> .*■«*'>'*-  ,2.17) 

where  we  choose  q  to  be  the  normal  coordinates  of  X3  and  choose  p  independently  of  q  so  that  7(q,p) 
is  only  a  function  of  q.  The  gas  phase  LEPS  potential  for  X3  can  be  conveniently  written  as  a  function 
of  the  X3  valence  coordinates  rlt  r2,  and  6,  in  which  r,  and  r2  are  the  bond  lengths  between  the  center 
atom  and  the  outer  two  atoms  respectively,  and  0  £  0  Srt  is  the  X-X-X  bend  angle.  Introducing  these 
variables  allows  us  to  evaluate  the  Jacobian,  and  thus  write  the  integration  over  internal  position  coordi¬ 
nates  in  Eq.  (2.17)  as 

J  dq  7(q)  e_*l(,y*'T  =  8jtV  j  dr,  j  tir2  J  </0  7(r„r *6)  i'a  W 

=  871^  |  dr 

where  V  is  the  volume.  If  we  write  an  "entropy-corrected"  potential  as 

<X>W2.e)  =  O^Wz*)  -  k„T  ln[J (r,/2,e)//(rj,r2,0*)]  ,  (2.19) 

in  which  removes  the  dimensionality  in  In,  then  Eq.  (2.18)  becomes 

J  dq  7(q)  e'*'Wk,T  =  7(rf,r$,0‘)  8rt*V  |  dr,  j  dr2  j  d9  .  (2.20) 

Figure  2  shows  the  shape  of  the  potentials  C>'(r,,r2,0),  C/£PS(rl,r2,Q)  and  the  function 
-ktT  In[7(r,,r2,0)/7(rf,r2,0*)]  versus  angle  0  for  r,  =  r2  =  2.254  A,  r\  =  r,,  i\  -  r2,  0*  =  90°  (i.e., 
sin(0*)  =  1),  and  LEPS  potential  parameters  defining  the  20  kcal/mol  barrier  (see  Table  2).  For  these 
parameters  the  minimum  in  energy  of  d>'(0)  occurs  at  0  -  148.5°.  The  use  of  internal  coordinate  vari¬ 
ables  and  the  effective  potential,  d>',  has  given  us  some  insight  into  the  nature  of  the  most  important 
position  coordinates  at  the  saddle  point;  we  now  see,  as  alluded  to  above,  that  the  bent  configuration  of 
X3  is  actually  more  probable  than  the  linear  configuration  of  X3  even  though  the  latter  is  a  lower  energy 
configuration  in  the  LEPS  potential.  The  effective  potential  "knows"  about  the  entropy  effect  of  having 
more  bent  configuration  space  available  to  compensate  for  the  higher  activation  energy  of  the  bent  X3 
geometries. 

The  effective  potential,  d>'(r,,r2,0),  is  used  in  the  selection  of  the  initial  position  and  momentum 
coordinates  for  X3.  Normal  modes  are  determined  at  the  potential  energy  minimum  of  <b'(r,,r2,0)  which, 
as  can  be  seen  from  Fig.  2,  is  well  approximated  by  a  quadratic  potential  close  to  the  minimum  in  the 
saddle  point  region.  These  normal  modes  may  be  interpreted  geometrically  as  a  bending  mode,  a  sym¬ 
metric  stretch  mode,  and  an  unstable  asymmetric  stretch  mode  which  has  imaginary  frequency  and  is 
used  to  define  the  reaction  coordinate.  Eq.  (2.11)  is  evaluated  by  selecting  values  of  r,,  r2,  and  0  using 
a  Gauss-Hermite  quadrature  in  the  normal  modes  of  G>\  as  suggested  by  Eq.  (2.20).  The  asymmetric 
stretch  variable  is  constrained  to  be  zero,  thus  defining  the  transition  state  surface  upon  which  all  trajec¬ 
tories  start  (It  is  not  even  a  good  idea  to  evaluate  the  integral  Eq.  (2.18)  by  the  above  method  using  the 
unmodified  potential,  due  to  the  fact  that  the  LEPS  potential  as  a  function  of  0  is  poorly 

represented  by  a  quadratic  potential  for  which  the  normal  mode  analysis  is  applicable.) 


|  dr2  |  dQ  raisin©  e~**£K(r'''i-6*i»T  ( 
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A s  previously  noted,  the  integration  over  the  momentum  variables,  p  in  Eq.  (2.17),  is  also  per¬ 
formed  using  Gauss-Hermite  quadrature  in  the  normal  modes,  and  the  momentum  quadrature  in  the 
asymmetric  stretch  is  chosen  to  be  even  to  ensure  that  there  will  always  be  some  initial  non-zero  velo¬ 
city  along  the  reaction  coordinate.  (Choosing  odd  quadrature  would  give  at  least  one  initial  condition 
for  which  the  asymmetric  stretch  momentum  is  zero).  For  every  selection  of  the  initial  conditions  of  X3, 
the  remaining  slow  position  coordinates  of  the  solvent  are  chosen  by  constraining  the  X3  geometry  at 
the  minimum  in  ®'(n/ 2,8)  while  allowing  the  solvent  to  equilibrate  using  molecular  dynamics.  This 
ensures  that  a  significantly  different  solvent  configuration  is  chosen  for  each  trajectory  in  the  ensemble. 
(Using  MD  in  this  manner  to  choose  initial  position  and  momentum  coordinates  for  the  solvent  requires 
considerable  computational  effort)  Finally,  the  momentum  variables  of  X3,  P,  are  selected  by  a  Gaus¬ 
sian  random  number  generator  in  a  Boltzmann  distribution  at  temperature  298  K. 

3.  Trajectories 

Trajectories  are  produced,  using  the  full  Hamiltonian,  by  propagating  both  forward  and  backward 
in  time  from  time  zero  at  the  initial  conditions  on  the  gas  phase  transition  state  surface.  For  trajectories 
which  turn  out  to  be  reactive,  propagating  in  one  direction  gives  the  dynamics  from  the  transition  state 
to  products  while  propagating  in  the  other  direction  gives  the  previous  time  history  from  reactants  to  the 
transition  state.  For  unreactive  trajectories,  both  forward  and  backward  propagation  leave  the  same 
atoms  bonded  together.  One  way  to  double  the  statistics  in  the  evaluation  of  Eq.  (2.11)  without  having 
to  select  twice  as  many  solvent  configurations  is  to  run  two  trajectories  for  each  solvent  configuration, 
one  trajectory  with  initially  positive  velocity  and  the  other  with  initially  negative  velocity  in  the  asym¬ 
metric  stretch. 

Finally,  we  have  verified  for  the  conditions  presented  here  that  there  is  no  recrossing  of  this  sur¬ 
face  in  the  gas  phase,  so  that  TST  is  exact  there. 

IQ.  X3  in  Argon 

In  this  section  we  present  our  results  for  the  special  case  of  the  hypothetical  reaction  of  X  +  X2  in 
liquid  argon,  with  X  having  the  mass  of  chlorine  35.  The  LEPS  parameters  are  initially  chosen  to 
describe  a  20  kcal/mol  barrier  height,  a  barrier  frequency  of  419  cm-1  and  the  isolated  molecule  Cl2 
Morse  potential36  (see  Table  2).  The  Lennard-Jones  parameters  (see  Table  1)  characterize  a  weakly- 
interacting  solvent  which  consists  of  100  argon  atoms  at  a  density  of  1400  kg  m-3  (per3  =  0.83) 
designed  to  approximate  the  density  of  argon  when  cool  enough  to  liquefy  at  a  pressure  of  1  atmo¬ 
sphere.  (The  actual  molecular  dynamics  are  carried  out  at  298  K.)  Results  for  this  system  described  in 
Sec.  A  are  then  analyzed  in  Sec.  B,  while  the  effects  of  lowering  the  barrier  frequency  are  considered 
in  Sec.  C. 

A.  Results 

According  to  simple  TST,  every  trajectory  which  arrives  at  the  transition  barrier  goes  directly  on 
to  form  products.  We  examine  this  assumption  for  the  20  kcal/mol  X3  in  argon  system  and  observe  that 
there  are  no  solvent  induced  barrier  recrossings  in  the  126  trajectories  sampled;  TST  is  exact.  The  rea¬ 
sons  for  this  are  described  below  in  Sec.  B. 

There  is,  however,  a  considerable  effect  of  the  solvent  on  the  details  of  the  reaction  trajectories  as 
they  approach  to,  and  recede  from  the  transition  barrier.  Figure  3  illustrates  how  the  solvent  can  perturb 
a  typical  reactive  trajectory  by  comparing  the  trajectory  on  the  gas  phase  potential  surface  and  in  solu¬ 
tion.  The  solvent  is  seen  to  change  the  precise  course  of  the  trajectory  away  from  the  barrier  compared 
to  its  unperturbed  path  in  the  gas  phase.  This  process  can  be  understood  in  simple  terms  by  examining 
the  time-dependent  exchange  of  energy  between  the  various  active  degrees  of  freedom  of  X3  and  the 
solvent  Using  Eq.  (2.11),  the  average  energy  of  X3  as  a  function  of  time  during  the  reaction  is  com¬ 
puted  for  126  trajectories.  After  the  center  of  mass  motion  of  X3  is  removed,  the  average  total  energy  of 
X3  is  partitioned  into  the  modes  of  X2  vibration,  rotation  and  translation  and  the  translational  energy  of 
the  X  atom.  This  partitioning  is  not  uniquely  definable  in  the  sense  that  other  definitions  are  possible, 
reflecting  the  lack  of  chemical  precision  in  defining  reactants  and  products  in  the  transition  hairier 
region  because  of  the  strong  X3  interaction  there.  We  choose  to  conveniently  define  the  two  closest  X 
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atoms  as  being  'diatomic"  X2  and  the  other  atom  as  being  "free".  When  one  atom  becomes  well 
separated  from  the  other  two  this  definition  becomes  chemically  meaningful  and  the  LEPS  potential 
reduces  to  just  the  bound  X2  Morse  potential.  The  total  "vibrational"  energy  is  then  taken  to  be  the 
sum  of  die  vibrational  kinetic  energy  of  the  closest  two  atoms  and  the  total  LEPS  potential  energy  of 
all  three  atoms.  (There  is,  of  course,  exchange  of  vibrational  energy  with  the  solvent)  Note  that  the 
zero  of  the  LEPS  potential  function  is  taken  to  be  the  potential  energy  of  an  isolated  X2  molecule  at  its 
equilibrium  geometry. 

Figure  4  shows  the  average  over  126  trajectories,  illustrating  the  time-dependent  partitioning  of 
energy  for  X3  in  argon  solvent  (solid  line)  at  T  «  298  K,  and  compares  this  with  the  average  X3  pardon¬ 
ing  as  a  function  of  time  in  the  gas  phase  (filled  circles).  Time  zero  is  at  the  transition  barrier.  The  plot 
is  symmetrical  about  time  zero  as  a  consequence  of  the  symmetry  in  the  potential  function.  It  is 
observed  that  there  is  a  short  period  of  time  before  and  after  the  transition  barrier,  between  approxi¬ 
mately  -0.05  ps  to  0.05  ps,  for  which  the  partitioning  of  energy  of  the  X3  in  solution  is  just  like  that  in 
the  gas  phase,  in  accord  with  the  representative  trajectory  shown  in  Fig.  3.  At  longer  times  away  from 
the  transition  barrier,  relaxation  of  the  X2  translational  and  rotational  energy  and  the  X  atom  transla¬ 
tional  energy  in  the  solvent  reduces  the  average  total  energy  of  the  X3  system  by  approximately  17 
kcal/mol  relative  to  the  initial  total  energy  at  the  transition  barrier.  No  such  decay  can  occur  in  the  gas 
phase  where  the  total  X3  energy  must  of  course  remain  constant.  The  X  and  X2  excess  translational 
energy  increases  sharply  as  one  falls  from  the  transition  barrier,  reaching  a  plateau  in  the  gas  phase  but 
rapidly  decaying  in  solution  as  energy  is  transferred  to  the  solvent.  These  processes  have  an  approxi¬ 
mately  02  ps  time  scale.  Similar  effects  can  be  seen  in  the  X2  rotational  energy  with  a  time  scale  of 
about  0.2  ps.  This  rotational  energy  is  small,  originating43  from  a)  the  repulsion  of  the  X  from  the  X2 
for  bent  original  configurations,  b)  bend  vibrational  motion  at  the  barrier,  and  c)  the  original  angular 
momentum  of  the  X3  at  the  transition  barrier. 

The  energy  in  the  separated  diatomic  X2  vibration  decays  in  approximately  0.25  ps  to  about  3  Has, 
where  to  is  the  ground  state  angular  vibration  frequency,  560  cm-1,  of  X2  =  CI2,  thus  making  the  pro¬ 
duct  diatomic  vibrationally  "hot"  for  a  comparatively  long  period  even  in  solution.  This  feature  is 
shown  by  Fig.  4  to  be  identical  in  both  phases  and  thus  is  governed  exclusively  by  the  LEPS  potential 
of  the  reaction  system.  One  expects  there  to  be  residual  vibrational  energy  in  part  because  the  X-X 
bonds  at  the  transition  barrier  are  slightly  longer  than  the  equilibrium  X2  bond  length.  This  excess 
vibrational  energy  does  not  relax  significantly  into  the  argon  solvent  in  a  1.0  ps  interval  (Fig.  4),  but 
rather  over  a  time  period  measuring  hundreds  of  picoseconds.  This  was  verified  by  integrating  one  typ¬ 
ical  trajectory  out  to  250  ps  and  observing  that  the  vibrational  energy  decays  over  this  time  period  to 
about  half  of  its  "initial"  value  of  3  Has. 

For  this  system,  the  reaction  outcome  is  determined  on  a  very  short  time  scale.  The  characteristic 
time  scale  for  the  reaction  fate  to  be  decided  is  less  than  0.1  ps  and  the  reaction  per  se  is  essentially 
finished  after  the  rotational  and  translational  energy  is  dissipated  to  the  solvent.  We  explore  this  aspect 
further  by  considering  the  space  and  time  dimensions  within  which  the  reaction  probability  is  sensitive 
to  solvent  momentum  perturbation  using  a  method  suggested  by  Andersen.29  An  ensemble  of  eighty 
reactive  trajectories  with  the  20  kcal/mol  barrier  in  argon  solvent  is  selected,  and  the  trajectories  are 
interrupted  at  some  time,  x,  before  crossing  the  barrier;  all  solvent  molecules  outside  a  sphere  of  a 
given  radius  measured  from  the  center  of  mass  of  X3  are  given  random  momenta  selected  from  a 
Boltzmann  distribution  while  keeping  all  position  coordinates  the  same.  Note  that  the  momenta  of  the 
reactive  trajectories  is  not  necessarily  Boltzmann.  The  ability  of  such  a  solvent  momentum  perturbed 
trajectory  to  cross  the  barrier  is  a  measure  of  the  solvent  momentum  "sphere  of  influence"  in  time  and 
space  upon  the  reaction.  Figure  5  shows  a  3D  plot  of  the  probability  of  successful  reaction  versus  both 
the  time,  x,  between  momentum  perturbation  and  original  trajectory  barrier  crossing,  and  the  radius  of 
the  spatial  sphere  outside  which  the  solvent  momentum  is  perturbed.  We  find,  as  one  might  intuitively 
expect,  that  for  a  constant  spatial  radius,  the  probability  of  reaction  decreases  for  increasing  time  while, 
for  constant  time,  the  probability  of  reaction  increases  for  increasing  radius.  If  all  of  the  solvent 
molecules’  momenta  are  perturbed  (i.e.  for  radius  set  equal  to  zero)  at  time  x,  then  it  is  seen  that  for  x 
>  0.14  ps  the  solvent  is  able  to  divert  all  trajectories  away  from  reaction.  Note  that  this  100%  diversion 
time  is  on  the  same  time  scale  as  the  fastest  energy  relaxation  process,  i.e.  translational  energy  relaxa¬ 
tion;  this  is  also  as  one  would  expect  on  the  basis  of  Fig.  4.  The  perturbation  in  solvent  momenta  must 
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evolve  by  translation  into  a  perturbation  in  solvent  positions  before  the  force  imparted  to  the  Xj  by  the 
solvent  is  altered,  and  thus  finally  the  motion  of  the  X3  over  the  barrier  can  be  deflected. 


B.  Discussion  and  Analysis 


1.  Recrossings 

We  observe  no  solvent  induced  recrossings  for  126  trajectories  with  the  X3  system  in  argon  with 
the  20  kcal/mol  barrier.  This  can  be  easily  understood  qualitatively  by  realizing  that  the  characteristic 
time  scale161*8  on  a  barrier  of  frequency  to*  -  419  cm-1  is  on  the  order  of  (to*)-1  -  0.01  ps.  On  this 
sharp  energy  barrier,  die  fate  of  the  trajectory  is  decided  very  quickly  (-  0.01  ps)  and  there  is  essen¬ 
tially  no  time  for  the  solvent  to  influence  the  outcome  of  the  reaction  through  any  effective  collisions 
with  the  X3.  This  idea  of  the  solvent  motion  and  its  effect  on  the  reaction  dynamics  can  be  expressed 
more  quantitatively  by  relating  the  time  dependent  solvent  friction  acting  on  the  reaction  coordinate16'17 
to  the  transmission  coefficient  k 


K  = 


k m 


(3.1) 


The  time  dependent  friction,  £(t),  is  given  by  the  time  correlation  function  of  the  fluctuating  forces  on 
the  reaction  coordinate 


C(0  =  0/H)  (  FF(i)  )  , 


(3.2) 


in  which  F  is  the  force  acting  on  the  reaction  coordinate  with  effective  mass  p  and  p  =  ( kBT)~ *.  This 
friction,  plotted  in  Fig  6,  is  computed  by  freezing  the  X3  coordinates  with  a  bent  X-X-X  geometry 
corresponding  to  the  minimum  in  potential  energy  of  the  effective  potential,  d>'(r,,r2,e)  (see  Sec.  II  and 
Fig.  2),  on  the  transition  barrier,  and  measuring  the  time  correlation  function  of  the  fluctuating  solvent 
forces  resolved  onto  the  reaction  coordinate.17  The  time  dependent  friction  is  also  computed  for  die  X3 
species  in  the  linear  geometry  on  the  transition  barrier  and  is  illustrated  in  Fig.  7.  This  friction  is  simi¬ 
lar  to  that  for  the  bent  geometry,  which  demonstrates  that  the  bend  angle  of  X-X-X  at  the  transition 
state  does  not  significantly  affect  the  nature  of  the  time  dependent  friction.  Because  of  the  very  fast 
passage  over  the  barrier,  the  effective  friction  acting  on  the  reaction  coordinate  is  not  the  full  time 
integrated  friction  £  given  by  Eq.  (1.2),  but  is  rather  closer  to  the  "instantaneous  friction"  felt  by  the 
reaction  system.  A  convenient  approximate  measure  of  this  is  the  value  of  the  friction  at  time  zero 


C(r  =  0)  =  (P/u)(F2)  m<a\. 


(3.3) 


which  we  describe  in  terms  of  an  initial  frictional  frequency,  ©5.  In  contrast,  the  full  friction,  £,  would 
be  important  for  long  time  scale  (e.g„  diffusion  controlled)  processes.  In  Fig.  6  we  see  that  the  solvent 
friction  for  argon  acting  on  the  bent  X3  system  falls  off  in  an  approximately  "Gaussian"  manner  with  a 
characteristic  time  t  -  0.19  ps  from  its  initially  weak  (©;  =  33  cm-1)  value;  in  particular,  ©t  is  small 
compared  to  the  barrier  frequency,  ©*,  value  of  419  cm-1.  Since  ©j  /  ©t  «  1  we  would  expect,16"18  as 
we  indeed  observe  in  the  dynamics,  dial  there  are  no  solvent  induced  recrossings  (in  126  trajectories), 
that  k  is  very  close  to  one,  and  thus  that  TST  is  valid.  This  result  will  also  hold  for  the  linear 
geometry  for  which  the  initial  value  of  the  friction  and  the  short  time  "Gaussian”  decay  is  the  same  as 
for  the  friction  of  the  bent  geometry  (see  Fig.  7).  (However,  the  time  decay  of  the  friction  is  slightly 
longer  for  the  linear  geometry  than  for  the  bent  geometry.  Thus,  for  the  short  time  scale  barrier  cross¬ 
ing,  the  linear  and  bent  X3  system  experiences  the  same  friction,  but  for  longer  time  scale  processes,  the 
linear  system  would  have  slighdy  more  friction  than  the  bent  system.) 


The  time  dependent  friction  on  the  reaction  coordinate  is  explicitly  incorporated  into  the  general¬ 
ized  Langevin  equation  theory  of  Grote  and  Hynes16"18  for  rate  processes.  Here,  k  is  calculated  as 


K—  X,  /  £0*  , 


(3.4) 


in  which  X,  is  termed  the  reactive  frequency  and  is  a  measure  of  the  solvent  response  to  the  reactive 
motion  across  a  barrier  of  frequency  ©*.  It  is  determined  from  the  self-consistent  equation 


K  = 


(3.5) 
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which  involves  the  hairier  frequency,  ©*,  and  die  Laplace  transform  frequency  component  of  the  time 
dependent  friction 

t(K)  =  |  dt  «"V  ;</)  (3.6) 

at  the  reactive  frequency  V  For  very  small  friction  Eqs.  (3.4)  and  (3.3)  give  X,  =  to*  and  k  =  1.  Also, 
for  sharp  barriers,  it  can  be  shown16*18  that  X,  =  to*  and,  again,  k  =  1.  If  we  consider  stronger  friction 
or  more  rounded  barriers,  then  it  is  expected  that  k  will  decrease  and  that  there  will  be  barrier  recross¬ 
ings.  With  a  Gaussian  approximation  to  the  decay  time  behavior  of  £(r)  (Fig.  6), 

C(r)  =  <&l  e*"*  (3.7) 

and  use  of  Eqs.  (3.4),  (3.3)  and  (3.6),  k  can  be  simply  determined  as  a  function  of  barrier  frequency  % 
for  various  values  of  (  and  the  results  are  shown  in  Fig.  8.  The  calculation  for  argon  with  b)»  -  419 
cm'1,  toj  =  35  cm*1,  and  the  relaxation  time  t  -  0.19  ps,  given  by  fitting  Eq.  (3.7)  to  the  data  in  Fig.  6, 
gives  k  m  1.0,  in  agreement  with  the  results  from  the  molecular  dynamics  computation. 

These  results  emphasize  that  it  is  the  solvent  response  at  the  reactive  frequency  which  determines 
how  effective  the  solvent  is  in  hindering  the  reaction  progress  over  the  barrier.  The  friction  felt  by  the 
reaction  system  during  the  time  period  (co*)-1  is  likely  much  smaller  than  the  long  time  friction  £.  Thus, 
it  is  not  surprising  that  the  20  kcal/mol  barrier  reaction  of  X3  in  argon  solvent,  which  spends  little  time 
on  the  barrier,  exhibits  very  few  recrossings  and  k  -  1.  Fig.  8  also  shows  that  for  smaller  values  of  ©t, 
k  is  predicted  to  eventually  decrease,  and  simple  TST  to  ultimately  fail.  We  explore  the  effect  of 
lowering  the  barrieT  frequency,  to*  in  Sec.  C  and  will  see  that,  for  lower  barrier  frequencies,  the  effect 
of  the  solvent  friction  increases  and  some  recrossings  are  observed. 


2.  Energy  Relaxation 

The  observed  basic  patterns  of  energy  relaxation  in  Fig.  4  are  what  we  intuitively  expect, 
Tntvtt  =  *Kor  «  * va,  and  can  largely  be  understood  in  terms  of  existing  models  of  energy  and  phase 
decay.44*47  Vibrational  energy  relaxation  by  the  solvent  is  much  slower  than  translational  and  rotational 
energy  decay  and  this  is  understood  by  considering  that  die  effective  friction  at  the  oscillator  frequency, 

Li  «*mc)  «  /  dt  e*osc  ‘  Ld(t)  ,  (3.8) 


estimated  by  a  Landau-Teller  type  expression  is  very  small  due  to  the  high  X2  oscillator  frequency  of 
approximately  560  cm*1.  The  "<T  subscript  emphasizes  that  the  time  dependent  friction  here  refers  to 
the  correlation  of  the  solvent  forces  on  the  diatomic  fragment. 

Since  the  remaining  3  Kto  vibrational  energy  in  the  nascent  product  X2  decays  only  slowly  and 
the  corresponding  amount  in  the  nascent  reactant  X2  is  only  acquired  slowly,  a  natural  question  arises: 
Is  vibrational  activation  up  to  £v  =  3  71©  of  the  reactants  the  rate  limiting  step?  If  so  then  the  rate  con¬ 
stant,  k,  will  depend  upon  the  dynamics  of  vibrational  activation  by  the  solvent  The  answer  is  no.  This 
can  be  explained  with  the  aid  of  the  Stable  States  Picture  of  reactions.10'18’48  For  a  symmetric  reaction 
the  influence  of  the  vibrational  activation  up  to  energy  Ev  is 


k  = 


k777 

1  +  TJ^Iky  ’ 


(3.9) 


in  which  the  influence  of  the  solvent  on  the  barrier  passage  is  justifiably  ignored  and  kv  is  the  vibra¬ 
tional  activation  rate  constant  up  to  Ev.  The  structure  of  kv  is 

kv  =  Zve'Uv,  (3.10) 


where  4  refers  to  the  details  of  the  vibrational  activation  process,  e  g.,  Landau-Teller  type  transitions. 
Clearly  vibrational  activation  dynamics  are  only  rate  limiting  if  they  represent  the  slow  step:  kv  «  kP7. 
But  since  £y  is  approximately  5  kcal/mol  while  the  activation  energy  is  20  kcal/mol,  we  instead  have 
knT  «c  e  &  «  ky.  Thus  the  barrier  passage  is  rate  limiting,  k  -» IP7  and  the  dynamics  of  vibrational 
energy  transfer  to  and  from  the  solvent  play  no  role  in  the  reaction  rate.  (The  relative  unimportance  of 
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vibrational  energy  transfer  in  atom  transfer  as  opposed  to  isomerizations  has  been  discussed10  else¬ 
where).  Not  only  must  the  X2  be  properly  vibrationally  excited,  but  also  the  liquid  must  be  in  a  part  of 
phase  space  where  it  can  transladonally  (and  rotationally)  thrust  the  X  atom  and  the  X2  molecule 
together  with  sufficient  energy  to  reach  the  barrier  top.  The  TST  rate  constant  reflects  the  equilibrium 
probability  of  all  of  these  events  occurring  simultaneously. 

Finally,  it  is  conceivable  that  the  approximately  17  kcal/mol  of  energy  rapidly  transferred  to  the 
solvent  in  the  reaction  could  have  resulted  in  a  local  "heating"  effect,  in  which  transladonally  energized 
solvent  atoms  quickly  reactivate  the  incipient  products  to  recross  the  barrier.  This  does  not  occur; 
instead  the  excess  energy  is  rapidly  transferred  to  the  other  solvent  atoms  outward  from  the  site  of  the 
reaction. 

3.  Simple  Model  for  Reaction 

A  simple  3-stage  picture  for  this  reaction,  as  shown  in  Fig.  9  emerges  in  light  of  the  above  obser¬ 
vations:  activation  of  reactants  in  stage  1,  barrier  crossing  in  stage  2,  and  deactivation  of  products  in 
stage  3.  In  the  first  stage  the  dynamical  role  of  the  solvent  is  to  furnish  activation  energy  to  the  reac¬ 
tants  and  properly  phase  their  motion  in  a  fashion  analogous  to  the  way  it  absorbs  the  excess  energy 
and  dephases  the  products  formed  in  stage  3.  There  are  three  time  scales  for  energy  exchange  with  the 
solvent  in  the  first  and  third  stages  corresponding  to  vibrational,  rotational  and  translational  energy 
exchange.  In  addition  to  solvent  induced  energy  activation  and  relaxation  (rt),  the  reaction  also  requires 
the  proper  phasing  among  relative  translational,  rotational  and  vibrational  motions  (r2).  The  second 
stage  includes  formation  from  die  energized  reacting  species  of  the  activated  complex  and  barrier  cross¬ 
ing.  Here,  for  to*  »  a> the  reaction  proceeds  over  the  barrier  top  largely  unimpeded  by  the  solvent, 
and  recrossings  are  not  observed.  In  Fig.  4  one  can  see  that  the  partitioning  in  the  liquid  and  gas  phases 
is  essentially  the  same  for  a  time  period  of  about  -0.03  ps  to  0.03  ps.  In  this  period  the  weakly- 
interacting  solvent  simply  does  not  have  enough  time  to  interfere  with  the  reaction  progress,  and  there 
is  essentially  unimpeded  gas  phase  dynamics  for  die  crossing  over  the  barrier.  This  statement  holds 
even  though,  as  shown  in  Fig.  3,  for  longer  times  before  and  after  the  transition  barrier  the  reaction  in 
the  liquid  phase  bears  little  resemblance  to  a  gas  phase  reaction.  As  soon  as  the  incipient  products  are 
formed  in  stage  3  they  rapidly  transfer  translational  and  rotational  energy  back  to  the  solvent,  and  are 
thus  energetically  prohibited  from  recrossing  the  barrier.  The  opposite  situation,  whereby  the  solvent 
transfers  additional  energy  to  the  products  after  they  are  formed,  is  not  observed  because,  as  discussed 
in  Sec.  m.B.2,  the  energy  of  the  products  is  high  compared  to  ktT  and  the  more  probable  event  is 
energy  deactivation.  Thus,  we  find  TST  to  be  an  excellent  description  for  this  reaction. 

We  now  examine  the  effect  on  the  dynamics  and  the  rate  constant  of  lowering  die  barrier  fre¬ 
quency. 

C.  Variation  of  Barrier 

In  Sec.  m.A  we  found  that  TST  is  appropriately  and  successfully  applied  in  the  region  -0.03  ps  to 
0.03  ps  during  which  time  the  barrier  is  crossed  and  the  fate  of  each  trajectory  is  decided.  We  find,  as 
predicted,16-18  that  TST  is  valid  for  the  case  of  a  sharp  barrier  and  a  weakly-interacting  solvent  But  if 
the  barrier  frequency  is  lowered  sufficiently,  then  the  characteristic  time,  (to*,)-1 ,  on  the  barrier  may 
become  sufficiently  long  that  the  solvent  has  time  to  induce  recrossings  (see  Fig.  8);  then  it  is  predicted 
that  TST  should  fail  to  some  extent  (k  <  1) 

Since  we  are  concerned  with  the  solvent’s  ability  to  disrupt  motion  across  the  barrier  and  the  sub¬ 
sequent  breakdown  of  TST,  we  will  consider  the  effect  of  changing  the  barrier  frequency,  to*,  as 
motivated  above.  A  simple  modification  of  the  LEPS  potential  parameters  iD„  i  -  1,2,3  (see  Table  2), 
allows  the  simultaneous  adjustment  of  the  barrier  height  and  barrier  frequency,  while  preserving  the 
asymptotic  characteristics.  The  barrier  height  has  no  direct  bearing  on  the  dynamical  influence  of  the 
solvent  although  it  will  affect  the  probability  of  formation  of  the  activated  complex  (and  hence  the  rate 

constant)  through  the  Boltzmann  weighting  /**'**r  where  E*  is  the  barrier  height  All  other  things  being 
equal,  a  lower  barrier  height  will  result  in  a  lower  barrier  freqiK  icy.  We  select  two  additional  barrier 
heights  of  10  kcal/mol  and  3  kcal/mol  with  corresponding  barrier  frequencies  288  cm-1  and  168  cm-1, 
as  shown  in  Table  2. 


Trajectories  are  run  on  the  above  potential  energy  surfaces  in  argon  solvent.  While  no  recross¬ 
ings  occur  for  the  10  heal/ mol  barrier,  some  recrossings  are  observed  for  the  S  keal/mol  barrier.  A  bar¬ 
rier  crossing  is  counted  each  time  that  the  value  of  the  asymmetric  stretch  reaction  coordinate  changes 
sign  during  die  course  of  the  trajectory,  beginning  with  an  initial  crossing  toward  either  the  product  side 
or  the  reactant  side.  Crossing  statistics  for  ensembles  of  trajectories  are  presented  for  the  20 , 10  and  5 
keal/mol  barriers  in  Tables  3a,  3b  and  3c  respectively.  The  statistics  in  the  tables  are  determined  by 
computing  trajectories  initialized  at  the  transition  barrier  with  either  positive  (toward  products)  or  nega¬ 
tive  (toward  reactants)  momentum  in  die  asymmetric  stretch  reaction  coordinate.  The  table  shows  that 
some  recrossings  occur  for  the  lowest  barrier  reaction. 


Table  3  (a).  Barrier  Crossings  for  20  keal/mol  X3  Barrier  in  Argon 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

126 

126 

0 

Table  3  (b).  Barrier  Crossings  for  10  keal/mol  X3  Barrier  in  Argon 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

32 

32 

0 

Table  3  (c).  Barrier  Crossings  for  5  keal/mol  X3  Barrier  in  Argon 


Initial  and  Final  States  of  Trajectories 


Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

27 

26 

0 

2 

5 

0 

0 

4 

3 

0 

0 

2 

0 

Thus,  Table  3c  indicates  that  while  27  trajectories  with  initially  positive  momentum  direedy 
formed  products,  5  trajectories  recrossed  quickly  to  form  reactants.  Similar  features  are  displayed  for 
trajectories  with  initially  negative  momentum.  However,  the  crossing  patterns  in  the  table  report  only 
the  number  of  trajectories  for  a  particular  number  of  crossings  and  give  little  indication  of  the  relative 
probabilities  for  those  trajectories  determined  by  the  quadrature  discussed  in  Section  Il.C.  One  estimate 
of  the  importance  of  these  probabilities  is  to  compute  the  average  number  of  crossings  with  proper 
weightings  using  Eq.  (2.11)  for  only  those  trajectories  which  react.  For  the  20  keal/mol,  10  keal/mol 
and  S  keal/mol  barriers  the  average  number  of  crossings  is  1.00,  1.00  and  1.03  respectively.  A  better 
estimate  for  die  recrossing  influence  on  the  rate  constant  is  the  transmission  coefficient  k. 

The  effect  of  the  recrossings  is  to  make  TST  an  overestimate  of  the  rate  constant,  and  the  correc¬ 
tion  to  the  TST  rate  constant  is  given  by  K  S  1.  We  will  determine  the  correction  k,  rather  than  the  rate 
constant  k  itself,  as  the  latter  requires  a  knowledge  of  the  free  energy  of  activation  due  to  solvent 
effects,  i.e.  the  potential  of  mean  force,  which  is  a  computationally  intensive  calculation  using  MD.49'51 
Furthermore,  K  is  a  direct  measurement  of  the  dynamical  influence  of  the  solvent  in  affecting  the  rate 
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constant. 

1.  Expression  for  x 

According  to  the  Stable  States  Picture  of  reactions48  the  rate  constant,  k,  can  be  expressed  by  the 
correlation  function  formula 
•• 

*  =  |*07*(0>a.  (3.11) 

in  which  j  is  the  flux  counted  positive  from  reactants  across  the  transition  barrier  at  time  zero  and  j*(t) 
is  the  flux  counted  negative  toward  products  at  a  later  time  t.  The  asterisk  indicates  that  the  dynamics 
of  the  trajectory  must  be  followed  only  long  enough  that  either  a  stable  reactant  or  product  is  formed 
(this  typically  occurs  in  our  reaction  system  in  much  less  than  1.0  ps)  and  (  )K  denotes  an  average  over 
the  equilibrium  distribution,  normalized  by  the  reactant  partition  function  QK.  The  time  integral  adds  up 
all  the  reactive  flux  into  stable  products  so  that  the  more  reactive  trajectories  there  are,  the  larger  is  the 
correlation  and  the  larger  is  the  rate  constant  k.  The  expression  for  i  in  Eq.  (3.11)  is  equivalently  writ¬ 
ten  upon  integration  as10 

*  =  lim  (  j  eW01  )a  .  (3.12) 

where  the  limit  means  we  follow  trajectories  only  long  enough  in  time  to  determine  whether  they  form 
stable  products  or  reactants,  0[x(/)]  is  a  step  function  which  equals  one  on  the  product  side  and  zero  on 
the  reactant  side  and  x(t)  is  the  asymmetric  stretch  reaction  coordinate  as  a  function  of  time  such  that  x 
is  zero  at  the  transition  barrier  for  r  «  0.  We  now  describe  various  calculational  forms  of  Eq.  (3.12) 
which  are  convenient  for  calculating  the  transmission  coefficient  K  and  closely  related  to  the  observed 
crossing  patterns  displayed  in  Table  3. 

In  simple  TST,  the  rate  constant  can  be  calculated  by  assuming  that  every  trajectory  at  the  transi¬ 
tion  barrier  which  has  initial  positive  flux  toward  products  will  always  end  up  as  products,  while  such 
trajectories  with  initial  negative  flux  toward  reactants  will  always  go  to  reactants.  Only  those  trajectories 
which  form  products  contribute  to  the  simple  TST  rate  constant,  and  there  are  never  any  recrossings. 
Thus,  in  calculating  t737,  we  need  to  consider  only  those  trajectories  which  have  initially  positive  flux, 
jM  and  there  is  no  need  to  follow  dynamics  for  any  time  past  zero  because  the  fate  of  all  trajectories  is 
decided  by  sign  of  their  initial  momenta.  The  flux  j+  is  given  by  the  velocity  across  the  transition  bar¬ 
rier,  (p/p)  5(x),  where  p  is  momentum,  p  is  the  reduced  mass  along  the  reaction  coordinate  of  the  X3 
species,  8(x)  is  the  delta  function  of  position  along  the  reaction  coordinate  with  x  -  0  at  the  barrier,  and 
the  TST  rate  constant  is4*6*20,21,39'52 

A”*- -<>♦>* 

=  (  (p/p)  8(x)  0[x(r)]  )„ 

=  (  (p/p)e(p))|,  (3.13) 

in  which  $  indicates  that  we  need  to  consider  only  the  initial  conditions  at  the  transition  barrier,  and 
0 (p)  is  a  step  function  which  is  one  for  positive  momenta  at  the  transition  barrier  and  zero  for  negative 
momenta. 

In  the  presence  of  recrossings,  the  actual  rate  constant,  Eq.  (3.12),  is  conveniently  split  into  two 
parts  representing  contributions  from  the  trajectories  with  initially  positive  momentum  and  those  with 
initially  negative  momentum  as 

k  =  (j  eWOl  )*  =  (/♦  0[x(r)] ),  +  (/_  0[x(»)]  ),■*.  +  *_,  (3.14) 

in  which  j+  represents  the  initially  positive  flux  while  j.  is  the  initially  negative  flux  across  the  transition 
barrier,  both  determined  by  the  magnitude  and  direction  of  the  initial  velocities  along  the  asymmetric 
stretch  reaction  coordinate  at  the  transition  barrier.  We  have  three  methods  of  applying  Eq.  (3.14)  to 
the  calculation  of  the  transmission  coefficient  for  the  rate  constant  using  MD. 

Method  1  is  to  sample  trajectories  which  always  have  initially  positive  velocities  at  the  transition 
barrier  at  time  zero.  This  focuses,  as  does  TST,  on  those  trajectories  initially  headed  toward  products 
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from  reactants.  For  the  contribution  k+,  die  fate  of  die  trajectories  is  determined  by  computing  the 
dynamics  forward  in  time.  For  the  contribution  k_,  we  imagine  time  reversing  the  velocities  of  all  par¬ 
ticles  in  the  reaction  system.  This  time  reversal  converts  j.  to  -7+  in  the  k_  contribution  to  k  in  Eq. 
(3-14)  as  well  as  changing  the  sign  of  /  in  ( 7.  8[x(r)]  )*.  Thus  k.  is  given  by 

(  7-  eWOJ  eW-t)]  >,  -  (3.15) 

The  step  function  0[x(-/)]  indicates  that  the  trajectories  must  originate  from  the  products  side  in  die 
asymptotic  past  Thus,  the  total  rate  constant  is  given  in  terms  of  initial  positive  velocity  trajectories  by 

k-(h  ®W0]  >*-(;♦  ew-r)]  )„  ,  (3.16) 

and  the  transmission  coefficient  k  is 

k  ( h  ewo]  )*  -  <  A  ew-r)] ), 

- nr, - ■  (3,7) 

Figure  10  illustrates  die  three  predominant  types  of  crossings  observed  (Table  3):  (a)  is  a  direct 
successful  reactant  -*  product  transition  with  no  recrossing;  (b)  is  a  single  recrossing,  reactant  -»  pro¬ 
duct  ->  reactant,  after  the  transition  barrier  is  initially  crossed  in  the  forward  direction,  and  (c)  is  a  sin¬ 
gle  reciossing,  product  -»  reactant  -»  product  in  which  the  recrossing  occurs  prior  to  a  reactant  -»  pro¬ 
duct  crossing  at  t  «  0.  The  equilibrium  distribution  on  the  transition  barrier  is  determined  by  the  algo¬ 
rithm  for  sampling  the  initial  conditions  using  Eq.  (2.11),  in  which  the  unnormalized  probability 
Wi  =  *T  and  the  initial  velocities  along  the  asymmetric  stretch,  are  determined  for  each  tra¬ 

jectory  i  in  the  ensemble  average.  Thus,  for  N  trajectories  in  the  ensemble  average,  k  is  computed  as 

Z  w‘  I  v»  I  Qi 

K  =  ,  (3.18) 

Z  k 1 

where  +  indicates  all  trajectories  have  initially  positive  velocity  across  the  transition  barrier,  and  the 
value  of  Qi  depends  on  the  initial  and  final  states  of  the  trajectory  i. 


if  Reactant 
if  Reactant 
if  Product  - 


■»  Product 

4  Reactant  or  Product 
Reactant 


Product 


As  described  in  Sec.  n.C.2,  before  a  particular  trajectory  is  computed,  the  Xj  coordinates  are  con¬ 
strained  while  solvent  coordinates  are  randomized  by  integrating  the  solvent  dynamics  out  to  many 
picoseconds.  This  ensures  that  a  significantly  different  solvent  configuration  is  chosen  for  each  trajec¬ 
tory  in  the  ensemble.  Using  MD  in  this  manner  to  choose  initial  position  and  momentum  coordinates 
for  the  solvent  requires  considerable  computational  effoit  One  way  to  double  the  statistics  in  the  deter¬ 
mination  of  k  without  having  to  select  twice  as  many  solvent  configurations  is  to  run  two  trajectories 
for  each  solvent  configuration,  one  trajectory  with  initially  positive  velocity  and  the  other  with  initially 
negative  velocity  in  the  asymmetric  stretch.  This  is  method  2  for  calculating  k,  in  which  the  trajec¬ 
tories  are  sampled  with  both  positive  and  negative  initial  asymmetric  stretch  velocities  and  time  reversal 
is  applied  to  all  trajectories  in  order  to  compute  the  dynamics  and  determine  the  initial  and  final  states. 
The  negative  flux  initial  conditions  are  chosen  from  the  positive  flux  initial  conditions  by  just  changing 
the  sign  of  the  velocities  of  the  Xj  while  keeping  all  position  coordinates  and  the  solvent’s  velocities 
unchanged.  Thus,  the  total  number  of  sampled  trajectories  is  twice  that  in  method  1  and  there  is  an 
equal  number  of  sampled  trajectories  from  both  the  positive  and  negative  flux  initial  conditions  with 
corresponding  equal  weights,  w„  and  equal  but  opposite  velocities,  v„  for  Xj.  Since  the  reaction  is  sym¬ 
metric  with  respect  to  reactants  and  products,  the  definition  of  reactants  and  products  is  a  matter  of  con¬ 
vention  for  a  given  trajectory.  In  method  2,  as  in  method  1,  we  want  to  sample  all  trajectories  which 
are  headed  initially  towards  reactants,  so  the  definition  of  reactants  and  products  will  depend  on 
whether  the  flux  is  positive  or  negative.  In  method  2,  all  positive  flux  initial  conditions  (also  used  in 
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method  1)  an  defined  as  heading  toward  products,  while  all  of  the  negative  flux  initial  conditions  an 
also  counted  as  heading  toward  products  for  the  TST  prediction  (only  for  method  2).  This  gives  k  as 


'n  I-1  T*  n 

:=  5>.iv.i  £ w, iv, i a  -  l w, |V, i e,  , 

i* *  L*>  <>- 


where  ±  indicates  all  trajectories  with  either  positive  or  negative  flux,  +  indicates  trajectories  with  posi¬ 
tive  flux  only  and  -  indicates  trajectories  with  negative  flux  only.  Eq.  (3.20)  can  be  shown  to  follow 
directly  from  Eq.  (3.19). 

Finally,  in  method  3,  we  use  the  positive  and  negative  flux  distributions  as  described  in  method  2 
but  keep  the  definition  of  nactants  and  products  the  same  for  both  the  positive  and  negative  flux  initial 
conditions.  This  convention  corresponds  to  the  definition  of  reactants  and  products  in  Table  3.  Thus,  all 
positive  flux  initial  conditions  are  headed  toward  products,  and  all  negative  flux  initial  conditions  are 
headed  toward  reactants.  Only  knowledge  of  the  final  outcome  of  the  reaction  and  the  magnitude  and 
direction  of  the  initial  flux  is  required  to  compute  the  transmission  coefficient  Trajectories  which  end 
up  as  products  contribute  either  positively  or  negatively  to  the  rate  constant  according  to  whether  their 
initial  flux  is  positive  or  negative.  This  leads  to  the  expression 


I  N  *  ’  N  N 

E  W,  |v,  I  £  W,  |v,  I  Q’,  -  |v,  |  Q'i  , 

•>  l*>+  <•- 

{1  if  Product 

0  if  Reactant  . 


J  1  if  Product 

0  if  Reactant  .  (3-22) 

This  last  formula  is  the  direct  implementation  of  Eq.  (3.14)  when  divided  by  k737.  It  is  also 
related  to  that  developed  by  Chandler53  and  employed  in  isomerization  studies.54  It  differs  in  that  tra¬ 
jectories  here  are  simply  followed  until  stable  reactants  or  products  are  formed,  rather  than  until  certain 
correlations  plateau  in  time.53'54 

All  three  methods  will  give  the  same  answer  in  the  limit  of  an  infinite  number  of  sampled  trajec¬ 
tories,  but  methods  2  and  3  are  ways  to  maximize  the  available  computer  time  by  doubling  the  number 
of  sampled  trajectories  without  having  to  recompute  a  new  set  of  initial  solvent  configurations. 
Methods  2  and  3  require  only  a  small  additional  amount  of  computer  time  compared  to  method  1. 

2.  Results  for  K 

For  the  S  kcal/mol  barrier  with  a  barrier  frequency  of  168  cm-1,  k  is  computed  from  an  ensemble 
of  64  MD  trajectories  and  using  methods  1,  2  and  3  as  0.91,  0.92  and  0.91  respectively.  Using  Eq. 
(3.4)  and  a  Gaussian  approximation  to  the  time  dependent  friction  for  argon  given  by  Eq.  (3.7)  with  x  - 
0.19  ps,  k  is  determined  from  Grote-Hynes  theory  as  0.98  for  the  5  kcal/mol  barrier.  For  the  10 
kcal/mol  barrier  (and  the  20  kcal/mol  barrier),  both  MD  and  Grote-Hynes  theory  predict  k  to  be  one. 
These  calculations  are  compared  in  Table  6.  The  agreement  in  the  full  MD  and  time  dependent  friction 
methods  of  calculating  k  is  excellent,  and  both  methods  give  the  lowest  value  of  k  for  the  lowest  bar¬ 
rier  frequency. 

It  is  clearly  very  interesting  to  go  to  lower  barriers  with  lower  a>»  values,  and  we  have  examined 
a  2.6  kcal/mol  barrier  which  has  a  barrier  frequency  of  cu*  -  20  cm-1.  However,  for  this  barrier  we  find 
multiple  recrossings  in  the  gas  phase.  These  apparently  arise  from  intramolecular  energy  flow  between 
the  reaction  coordinate  and  coupled  non-reactive  modes.55  The  longer  time  spent  on  the  barrier,  (to*,)-1 
-  0.3  ps,  may  allow  such  coupling  effects  to  develop.  These  isolated  molec  ile  dynamics  must  be  under¬ 
stood  before  the  solvent  influence  on  the  reaction  can  be  properly  characterized,  and  this  topic  is  left  for 
future  research.  Nonetheless  it  is  worth  stressing  this  observation  as  a  cautionary  example  of  the  possi¬ 
ble  unphysical  character  of  low-banier,  one-dimensional  reaction  studies  which  ignore  other  internal 
degrees  of  freedom. 


m 
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3.  EaertJ  Relaxation 

The  eneigy  decay  patterns  for  X5  in  argon  for  the  10  kcal/mol  and  S  kcal/mol  barriers  are  similar 
to  those  for  the  20  kcal/mol  barrier  in  Fig.  4  but  are  scaled  down  in  magnitude  in  proportion  to  the  bar¬ 
rier  heights.  The  eneigy  decay  patterns  are  presented  for  the  5  kcal/mol  barrier  in  Fig.  11.  An  interest¬ 
ing  new  feature  is  the  transient  rise  in  solution  of  the  potential  energy  plus  diatomic  vibrational  kinetic 
energy  of  X2  above  the  gas  phase  value.  We  have  traced  this  to  a  certain  "cage"  effect  in  which  the 
high  density  argon  solvent  temporarily  holds  the  diatomic  fragment  in  a  region  of  finite  X3  potential 
energy. 

IV.  Variation  of  Solvent 

In  Sec.  ID  above,  calculations  of  energy  decay,  solvent  momentum  sphere  of  influence,  recross¬ 
ings,  and  k  were  presented  for  the  reaction  X  +  X2  in  argon  solvent  In  this  section  we  briefly  explore 
the  corresponding  results  in  helium  and  xenon  solvents.  These  represent  "light"  and  "heavy"  molecule 
solvents  when  viewed  on  the  mass  scale  of  X  -  Cl. 

The  solvent  systems  consist  of  100  atoms  at  densities  of  147  kg  m-3  (po3  =  0.36)  and  3520  kg 
m~3  (per3  =  0.80)  for  helium  and  xenon  respectively,  designed  to  represent  the  density  of  these  solvents 
at  a  temperature  cool  enough  to  liquefy  at  a  pressure  of  1  atmosphere.  The  actual  molecular  dynamics 
are  computed  for  a  temperature  of  298  K. 

A.  Energy  Relaxation 

The  patterns  for  eneigy  relaxation  are  shown  in  Figs.  12  and  13  for  the  20  kcal/mol  eneigy  bar¬ 
rier  in  helium  and  xenon  respectively.  The  corresponding  results  for  argon  are  shown  in  Fig.  4.  The 
translational  and  rotational  energy  decay  times  are  faster  in  argon  and  xenon  solvents  and  slower  in 
helium  solvent  Helium,  with  a  small  mass,  should  be  iess  effective  in  exchanging  these  types  of 
eneigy  through  collisions  with  X  and  X2,  and  indeed  energy  decay  to  and  from  translational  and  rota¬ 
tional  motion  is  relatively  slow.  Xenon,  with  a  large  mass,  is  similar  to  argon  in  its  ability  to  dissipate 
the  excess  translational  and  rotational  energy.  These  patterns  are  not  unexpected.  Although  this  is  not 
shown  in  the  figure,  only  helium  begins  to  relax  vibrational  energy  in  less  than  10  ps.  This  is  in  rough 
accord  with  a  Landau-Teller  picture  in  which  high  frequency  force  components  arising  from  rapid 
helium  motions  are  important 

As  with  argon,  the  energy  decay  patterns  for  lower  banier  heights  are  similar  to  those  for  the  20 
kcal/mol  barrier  but  scaled  down  in  magnitude.  The  energy  relaxations  for  the  5  kcal/mol  barrier  in 
helium  and  xenon  solvents  are  illustrated  in  Figs.  14  and  15  respectively.  The  cage  effect  discussed  in 
Sec.  IV.C.3  for  the  transient  excess  of  Xj  potential  plus  diatomic  vibrational  kinetic  energy  for  the  5 
kcal/mol  banier  in  argon  appears  again  in  Fig.  15  for  the  high-density,  heavy  solvent  xenon,  but  is 
absent  in  Fig.  14  for  the  lower  density,  light  solvent  helium. 

Finally,  the  solvent  momentum  sphere-of-influence  for  the  20  kcal/mol  barrier  in  helium  solvent 
is  shown  in  Fig.  16.  The  momentum  effect  on  the  reaction  by  helium  solvent  is  lower  in  magnitude 
than  for  argon  in  Fig  5. 

B.  Recrossings 

Trajectories  are  analyzed  for  recrossing  patterns  and  the  results  for  the  20,  10  and  5  kcal/mol 
energy  barriers  in  helium  and  xenon  solvents  are  presented  in  tables  4a,  4b,  4c,  5a,  5b,  and  5c.  There 
are  no  recrossings  observed  for  the  20  kcal/mol  and  10  kcal/mol  barriers  in  any  of  the  three  solvents. 
Also  note  that  there  are  no  recrossings  observed  on  any  of  the  barriers  in  the  absence  of  solvent  (i.e.  in 
the  gas  phase).  The  average  number  of  crossings  is  computed  for  the  5  kcal/mol  barrier  in  helium  and 
xenon  as  1.00  and  1.01,  respectively.  This  average  includes  the  relative  weights  of  the  trajectories. 

In  Sec.  in  we  used  the  idea  of  the  time  dependent  solvent  friction  acting  on  the  reaction  coordi¬ 
nate  and  found  that,  for  the  20  kcal/mol  barrier,  the  initial  frictional  frequency  ci);  in  argon  solvent  is 
much  weaker  than  the  frequency  associated  with  the  forces  driving  the  reaction  down  the  barrier,  that 
recrossing  are  accordingly  rare,  and  that  the  transmission  coefficient  k  =  1.  For  the  low,  5  kcal/mol 
barrier  height  and  frequency,  some  recrossings  are  observed  and  tc  decreases  slightly.  The  two 
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Table  4  (a).  Barrier  Crossings  for  20  kcal/moi  X3  Barrier  in  Helium 


Number  of  Crossings 


1 


Initial  and  Final  States  of  Trajectories 


Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

0 

32 

32 

0 

Table  4  (b).  Barrier  Crossings  for  10  kcal/moi  X3  Barrier  in  Helium 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  •  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

32 

32 

0 

Table  4  (c).  Barrier  Crossings  for  5  kcal/moi  Xj  Barrier  in  Helium 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

30 

31 

0 

2 

1 

0 

0 

2 

additional  solvents  helium  and  xenon  have  a  time  dependent  friction  which  diffeis  from  argon,  as  illus¬ 
trated  in  Figs.  6  and  7.  Helium  has  a  weaker  friction  than  argon  with  initial  frequency  co^  =  15  cm-1, 
and  a  short  decay  time  T  of  about  0.06  ps,  while  the  xenon  friction  is  stronger  with  b>(  =  40  cm-1  and 
decay  time  x  of  0.36  ps.  We  observe  that  recrossings  in  these  solvents  (see  Tables  4  and  5)  for  .he  low 
barrier  of  5  kcal/moi  are  similar  to  what  is  seen  for  the  low  barrier  height  in  argon  (see  Table  3). 
There  are  however,  compared  to  argon,  fewer  recrossings  in  helium  but  more  recrossings  in  xenon.  This 
is  in  accord  with  the  calculation  that  there  is  less  short  time  solvent  friction  for  helium  and  greater  short 
time  friction  for  xenon  in  comparison  to  argon.  None  of  the  solvents  presents  a  great  enough  friction 
on  the  reaction  coordinate  to  induce  recrossings  for  barrier  heights  of  10  and  20  kcal/moi,  as  summar¬ 
ized  in  Table  6. 

Using  methods  1,  2  and  3  as  described  in  Sec.  Ill,  k  is  computed  for  64  trajectories  on  the  5 
kcal/moi  barrier  in  helium  as  0.96,  0.97  and  0.98,  respectively,  while  k  for  252  trajectories  on  the  5 
kcal/moi  reaction  in  xenon  is  0.91,  0.91,  and  0.92.  As  seen  here  and  in  Table  6,  k  decreases  for  lower 
barrier  heights  and  for  stronger  solvent  time  dependent  frictions.  This  is  compared  with  Grote-Hynes 
(GH)  theory,  Eqs.  (3.4)  -  (3.6),  for  computing  k,  when  we  use  a  Gaussian  approximation  to  the  time 
dependent  friction.  For  helium  with  a  5  kcal/moi  barrier,  barrier  frequency  to*  -  168  cm-1  and  T  «  0.06 
ps,  GH  theory  gives  k  as  0.99,  while  for  xenon,  with  t  -  0.36  ps,  k  is  predicted  to  be  0.95.  A  com¬ 
parison  of  k  computed  from  MD  and  GH  theory  is  presented  for  all  of  the  barrier  heights  and  solvents 
in  Table  6.  The  estimated  error  in  the  values  of  k  from  MD  for  the  5  kcal/moi  barrier  is  ±0.03.  The 
values  for  k  from  GH  theory  compare  within  0.05  of  the  values  determined  by  molecular  dynamics 
simulation  and  the  trend  for  k  to  decrease  with  decreasing  barrier  frequency  and  increasing  initial  sol¬ 
vent  frequency  are  the  same  for  both  GH  theory  and  MD. 

The  5  kcal/moi  reaction  in  xenon  is  an  interesting  one  with  which  to  calculate  the  predictions  of 
Kramers  Theory.19  Since  die  xenon  time  dependent  friction  is  of  significant  magnitude  and  fairly  long 
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Table  5  (a).  Barrier  Crossings  for  20  kcal/mol  X3  Barrier  in  Xenon 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  •  Reactant 

Product  -  Product 

1 

0 

126 

126 

0 

Table  5  (b).  Barrier  Crossings  for  10  kcal/mol  X3  Barrier  in  Xenon 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

32 

32 

0 

Table  5  (c).  Barrier  Crossings  for  5  kcal/mol  X3  Barrier  in  Xenon 

Number  of  Crossings 

Initial  and  Final  States  of  Trajectories 

Reactant  -  Reactant 

Reactant  -  Product 

Product  -  Reactant 

Product  -  Product 

1 

0 

72 

74 

0 

2 

19 

0 

0 

21 

3 

0 

31 

21 

0 

4 

7 

0 

0 

4 

5 

0 

1 

2 

0 

lifetime,  (Figs.  6  and  7),  the  friction  constant 


;  =  (P/H)J  dt  {  FF(/) )  ,  (4.1) 

is  considerable.  We  numerically  estimate  from  Figs.  6  and  7  that  £  =  105  cm"’  for  the  bent  X3 
configuration  and  £  =  198  cm-1  for  the  linear  configuration.  Then  the  Kramers  transmission 
coefficients19 

k  =  V 1  +  (£/2(i>t)2  -  (C/2o>*)  (4.2) 

are  0.74  and  0.59  for  the  bent  and  linear  configurations.  These  are  noticeably  less  than  the  MD  results 
and  the  Grote-Hynes  prediction,  emphasizing  again  the  importance  of  the  short  time  scale  solvent 
dynamics  for  the  reaction  problem. 

This  short  time  aspect  can  be  finally  stressed  by  completely  ignoring  the  time  dependence  of  the 
friction.  Then  £(X,)  in  Eq.  (3.5)  reduces  to  to|  /  X,  and  we  obtain  the  short  time,  nonadiabatic  prediction 
of  GH  theory16-1* 

1C  =  Vl  -(<■>;/  to»)2  (4.3) 

solely  in  terms  of  the  barrier  frequency  and  the  initial  time  dependent  friction  value.  Equation  (4.3) 
agrees  to  within  0.01  with  the  Gaussian  friction  GH  k  values  listed  in  Table  6  and  is  thus  in  very  good 
agreement  with  the  MD  k  values. 
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Table  6. 

Transmission  Coefficients  for  Various  Barriers  and  Solvents 

Barrier 

Barrier 

K 

K 

Average  Number  of 

Height 

(kcal/mol) 

Frequency 

(cm-1) 

Solvent 

(Molecular  Dynamics) 

(Grote-Hynes  Theory) 

Barrier  Crossings 

He 

1.00 

1.00 

1.00 

20 

419 

Ar 

1.00 

1.00 

1.00 

Xe 

1.00 

0.99 

1.00 

He 

1.00 

1.00 

1.00 

10 

288 

Ar 

1.00 

0.99 

1.00 

Xe 

1.00 

0.98 

1.00 

He 

0.97 

0.99 

1.00 

5 

168 

Ar 

0.91 

0.98 

1.03 

Xe 

0.91 

0.95 

1.01 

Finally,  the  three  solvents,  helium,  argon  and  xenon,  display  weak  coupling  with  the  reaction 
coordinate.  This  is  reflected  in  the  lack  of  recrossings  observed  for  the  higher  barriers  and  relatively  few 
recrossings  on  the  lowest  barrier.  A  more  strongly  interacting  solvent,  which  has  a  larger  initial  value  of 
the  friction,  ©j,  would  be  expected  to  decrease  k.  One  way  to  increase  the  solvent-solute  coupling 
strength  is  to  modify  the  value  of  the  X3-solvent  Lennard-Jones  parameter,  When  the  value  of  ey 
for  argon-X  interactions  is  multiplied  by  a  factor  of  128  we  find  that  the  dynamics  with  the  20  kcal/mol 
are  not  significantly  different  than  the  dynamics  for  the  unchanged  solvent  The  rate  of  energy  exchange 
with  this  modified  solvent  is  much  faster,  but  the  transition  barrier  dynamics  are  the  same;  there  are  no 
recrossings,  and  k  is  one.  If  the  solvent  interaction  parameters  are  modified  in  this  manner  the  solvent 
becomes  very  "sticky"  and  stays  roughly  locked  in  a  particular  configuration.  This  slows  down  the  sol¬ 
vent  motion  significantly,  and  barrier  passage  occurs  for  a  solvent  configuration  which  is  essentially 
fixed.  However,  even  for  this  major  increase  in  the  solvent-solute  coupling  strength,  the  dynamics  on 
the  barrier  are  still  roughly  approximated  by  gas  phase  barrier  dynamics.  Unfortunately,  the  "sticky" 
character  of  the  solvent  prevented  us  from  obtaining  reliably  equilibrated  transition  barrier  initial  condi¬ 
tions  sufficient  to  pursue  this  case  quantitatively. 


V.  Conclusions 

Molecular  dynamics  have  been  computed  for  the  hypothetical  reaction  of  X  +  X2  -*  X2  +  X  (the 
mass  of  X  set  to  chlorine)  in  three  different  liquid  density  inert  gas  solvents,  He,  Ar,  and  Xe.  We  find 
that  this  reaction  can  be  largely  understood  in  terms  of  a  simple  model  of  3  epochs:  a)  build  up  of  the 
reactant  energy  and  phase  relationships  needed  for  activation,  b)  barrier  crossing,  and  c)  decay  of 
energy  and  phase  of  products.  Epochs  a)  and  c)  can  largely  be  understood  in  terms  of  existing  models 
of  phase  and  energy  decay  while  b)  can  largely  be  understood  in  terms  of  gas  phase  A  +  BC  barrier 
crossing  dynamics,  even  for  the  lowest  barriers  we  have  examined. 

We  find  no  solvent  induced  recrossings  and  thus  no  deviation  from  TST  for  reactions  with  higher 
barriers  and  barrier  frequencies.  A  very  small  deviation  from  gas  phase  barrier  crossing  dynamics  and 
thus  a  slight  breakdown  of  TST  due  to  recrossings  of  the  barrier  is  found  upon  lowering  the  barrier 


frequency;  the  characteristic  time  spent  at  the  top  of  the  barrier  becomes  longer  and  somewhat  more 
comparable  to  the  time  scale  for  solvent  motion,  and  the  effect  of  friction  by  the  solvent  on  the  reaction 
coordinate  increases.  The  model  of  time  dependent  friction  acting  on  the  reaction  coordinate  is  shown 
to  be  useful  in  understanding  and  quantitatively  predicting  the  role  of  the  solvent  in  affecting  the 
dynamics  of  barrier  crossing  and  therefore  the  rate  constant. 

Three  methods  for  computing  the  transmission  coefficient  from  molecular  dynamics  have  been 
presented  and  the  computed  results  compared  to  the  predictions  of  the  Grote-Hynes  equations  in  which 
the  initial  time  dependent  friction  is  approximated  by  a  Gaussian.  The  molecular  dynamics  and  Grote- 
Hynes  theory  give  values  for  the  transmission  coefficient,  k,  which  compare  within  0.0 5  and  give  the 
same  trends  of  decreasing  k  with  increasing  solvent  friction  and  decreasing  barrier  frequency. 

We  have  discussed  die  very  low  reaction  barrier  (2.6  kcal/mol)  case  in  only  a  cursory  fashion 
here,  due  to  the  apparent  involvement  of  intramolecular  energy  flow  features  in  the  isolated  gas  phase 
reaction,  thereby  clouding  the  role  of  the  solvent  We  hope  to  analyze  this  interesting  case  in  detail  in 
the  future.  Other  interesting  possibilities  for  future  study  include  strongly  interacting  solvents  involving, 
for  example,  hydrogen  bonding  to  the  reaction  system.56*57 

We  have  studied  here  only  the  symmetric  reaction  case.  There  are  new  dynamical  aspects  present 
in  asymmetric  reactions;  for  example,  the  role  of  reactant  vibrational  excitation  in  accelerating  reac¬ 
tions.7  It  will  be  interesting  to  examine  whether  these  features  survive  in  solution.  Our  present  results 
are  already  highly  suggestive  in  this  regard. 

A  major  finding  of  this  work  is  that  TST  typically  provides  an  excellent  description  for  the  rates 
of  simple  atom  transfers  in  weakly  interacting  solvents.  This  means,  for  example,  that  one  can  with 
confidence  use  the  methods  of  modem  equilibrium  statistical  mechanics  to  compute  solution  rate  con¬ 
stants  using  TST  for  such  reactions.  In  this  connection  it  should  also  be  stressed  that  in  addition  to  the 
dynamical  effects  discussed  here,  within  the  TST  description  there  can  be  important  equilibrium  solvent 
effects  on  atom  transfer  free  energy  barriers  and  thus  on  rate  constants.58  The  present  work  suggests 
that  when  these  are  accounted  for,  one  might  reliably  predict,  for  example,  the  pressure  variation  of  a 
simple  atom  transfer  rate  from  gas  to  liquid  phase  densities  in  inert  rare  gas  solvents. 
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Figure  Captions 


Figure  1.  Potential  energy  LEPS  surface  contour  plot  as  a  function  of  the  AB  and  BC  bond  lengths,  r, 
and  r2,  for  a  linear  arrangement  of  X3  with  a  barrier  height  of  20  kcal/mol,  and  with  a  barrier  frequency 
of  to*  -  419  cm-1.  The  saddle  point  is  indicated  by  a  solid  circle.  The  contour  lines  are  labeled  in  units 
of  kcal/mol.  Note  that  the  zero  of  the  LEPS  potential  function  is  taken  to  be  the  potential  energy  of  an 
isolated  X2  molecule  at  its  equilibrium  geometry. 

Figure  2.  Effective  potential  d>'(8),  ^£PS(6)  and  the  function  -kBT  //i[/(8)/rfrf]  as  a  function  of  8  for 
the  20  kcal/mol  barrier  of  X3  in  which  rx  and  r2  are  held  constant  at  2.254  A. 


Figure  3.  Trajectory  plots  of  X3  on  the  20  kcal/mol  LEPS  potential  surface  for  the  reaction  in  both 
liquid  xenon  (dashed  line)  and  die  gas  phase  (solid  circles)  with  identical  initial  conditions  chosen  at 
time  zero  on  die  transition  state  surface,  demonstrating  the  solvent  effect  on  the  details  of  the  dynamics 
away  from  the  transition  barrier.  The  saddle  point  of  the  potential  is  indicated  by  a  single  solid  circle. 


Figure  4.  Plot  of  die  average  energy  versus  dme  for  the  reaction  X  +  X2-»X2  +  X  with  the  20 
kcal/mol  energy  barrier,  for  an  ensemble  of  126  reacting  trajectories.  Time  zero  is  when  the  X3  is  first 
released  at  the  saddle  point  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  argon 
solvent  while  the  dots  indicate  the  average  energy  for  reaction  in  the  absence  of  solvent  (gas  phase). 
The  translational  energy  arises  from  and  decays  to  the  solvent  in  0.2  ps,  rotational  in  about  0.2  ps  and, 
although  the  figure  does  not  show  this,  vibrational  energy  in  >  250  ps.  Notice  how  translational  and 
rotational  energies  are  dumped  into  ABC  potential  energy  to  climb  the  steep  barrier,  and  then  dumped 
back  into  translation  and  rotation.  The  vibrational  energy  of  the  X2  before  and  after  barrier  passage  is 
approximately  3  Hto,  where  to  is  the  560  cm-1  ground  state  vibrational  angular  frequency  of  X2.  Note 
that  the  average  energies  in  the  different  modes  for  the  gas  and  liquid  phase  reactions  are  approximately 
the  same  during  the  time  period  -0.05  ps  to  0.05  ps  during  which  the  gas  and  liquid  phase  reaction  tra¬ 
jectories  are  similar. 


Figure  5.  Solvent  momentum  sphere  of  influence  on  reaction  with  20  kcal/mol  barrier  in  argon  solvent 
For  a  chosen  ensemble  of  80  trajectories  which  react,  the  trajectories  are  carried  back  the  period  in  time 
before  the  barrier  crossing  given  on  the  horizontal  axis.  While  holding  all  coordinates  constant,  all  sol¬ 
vent  atoms  outside  a  sphere  of  radius  measured  from  the  center  of  mass  of  the  X3  atoms  are  random¬ 
ized  with  a  Boltzmann  distribution  in  velocities.  The  trajectories  are  then  allowed  to  resume.  The  pro¬ 
bability  on  the  vertical  scale  measures  the  fraction  of  trajectories  which  still  lead  to  reaction.  This  gives 
a  measure  of  the  time  and  space  scale  over  which  solvent  momenta  play  a  role  in  the  reaction. 


Figure  6.  Solvent  friction  on  the  reaction  coordinate  for  the  X3  in  a  bent  geometry  (defined  by  the 
minimum  in  the  effective  potential,  d>\  in  Fig.  2)  as  function  of  time  for  the  three  inert  solvents  helium, 
argon  and  xenon,  normalized  by  the  value  of  the  friction  at  time  zero.  The  initial  decay  in  the  time 
dependent  friction  is  approximated  by  a  Gaussian  Eq.  (3.7)  with  values  for  t  of  0.06  ps,  0.19  ps,  and 
0.36  ps  for  He,  Ar  and  Xe  respectively,  determined  by  numerically  calculating  the  curvature  of  the  time 
dependent  friction  at  I  -  0.  The  initial  frictional  frequencies,  CD;,  for  the  three  solvents  are  15  cm'1,  35 
cm"1  and  40  cm"1  for  He,  Ar  and  Xe  respectively. 
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Figure  7.  Solvent  friction  on  the  linear  geometry  reaction  coordinate  as  function  of  time  for  die  three 
inert  solvents  helium,  argon  and  xenon,  normalized  by  the  value  of  the  friction  at  time  zero.  The  initial 
frictional  frequencies  for  the  three  solvents  as  well  as  the  initial  decay  in  the  dependent  friction  are 
essentially  the  same  as  for  the  bent  case  in  Fig.  6.  The  time  dependent  friction  for  the  linear  geometry, 
however,  has  a  longer  time  decay  than  that  for  the  bent  geometry  shown  in  Fig.  6. 


Figure  8.  The  value  of  the  transmission  coefficient,  k,  from  Grote-Hynes  theory  as  a  function  of  barrier 
frequency,  to*  in  die  solvents  helium,  argon,  and  xenon.  The  initial  decay  in  the  the  time  dependent 
frictions  for  these  solvents  (see  Fig.  6)  can  be  approximated  by  Gaussians  Eq.  (3.7),  with  values  for  t  of 
0.06  ps,  0.19  ps,  and  0.36  ps  for  He,  Ar  and  Xe  respectively.  The  Gaussian  approximation  is  good  for 
barrier  frequencies  greater  than  SO  cm-1  where  the  time  spent  on  the  barrier  is  short,  but  poor  for  bar¬ 
rier  frequencies  less  than  50  cm-1  since  die  longer  time  friction  becomes  more  important  as  the  charac¬ 
teristic  time  spent  on  the  barrier  increases.  For  the  reaction  systems  described  in  this  work,  the  short 
time  Gaussian  approximation  is  adequate. 


Figure  9.  Illustration  of  the  simple  3-stage  picture  to  explain  the  essential  features  of  the  dynamics  of 
A  +  BC  reactions  in  rare  gas  solvents.  The  picture  consists  of  three  stages.  Stage  1  is  energy  and  phase 
arisal  in  reactants,  stage  2  is  approximately  gas  phase  barrier  crossing,  while  stage  3  is  energy  and 
phase  decay  of  products  in  the  solvent 


Figure  10.  Schematic  illustration  of  some  predominant  recrossing  patterns  observed  in  the  MD  simula¬ 
tion  of  A  +  BC  reaction  dynamics  on  low  barriers  in  solution.  The  dividing  line  represents  the  transition 
barrier  dividing  surface  between  reactants  on  the  left  and  products  on  the  right  while  die  arrows  indi¬ 
cate  the  direction  of  the  trajectories  for  initially  positive  momentum:  (a)  is  a  direct  successful  reactant 
-»  product  transition  with  no  recrossing;  in  TST,  all  initially  forward  trajectories  are  of  this  type;  (b)  is 
a  single  recrossing,  reactant  -»  product  ->  reactant  after  the  transition  barrier  is  crossed  in  the  forward 
direction,  and  (c)  is  a  single  recrossing,  product  -♦  reactant  ->  product  in  which  the  recrossing  occurs 
prior  to  a  reactant  — >  product  crossing. 


Figure  11.  Plot  of  the  average  energy  versus  time  for  the  reaction  X  +  X,  ->  X2  +  X  with  the  5 
kcal/mol  energy  barrier,  for  an  ensemble  of  64  trajectories.  Time  zero  is  when  the  Xj  is  first  released  at 
the  saddle  point  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  argon  solvent  while 
the  dots  indicate  the  average  energy  for  a  reaction  in  the  absence  of  solvent  (gas  phase).  The  transla¬ 
tional  and  rotational  energies  arise  from  and  decay  to  the  solvent  in  approximately  0.5  ps,  and,  although 
the  figure  does  not  show  this,  the  vibrational  energy  arises  and  decays  in  >  250  ps. 


Figure  12.  Plot  of  the  average  energy  versus  time  for  the  reaction  X  +  X2  — >  X2  +  X  with  the  20 
kcal/mol  energy  barrier,  for  an  ensemble  of  64  reacting  trajectories.  Time  zero  is  when  the  X3  is  first 
released  at  the  saddle  point.  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  helium 
solvent  while  the  dots  indicate  the  average  energy  for  a  reaction  in  the  absence  of  solvent  (gas  phase). 
The  translational  and  rotational  energy  arises  from  and  decays  to  the  solvent  over  a  time  scale  measur¬ 
ing  longer  than  1.0  ps,  and,  although  die  scale  of  the  figure  does  not  show  it,  vibrational  energy  starts  to 
decay  appreciably  in  less  than  10  ps. 


Figure  13.  Plot  of  the  average  energy  versus  time  for  the  reaction  X  +  X2  -»  X2  +  X  with  the  20 
kcal/mol  energy  barrier,  for  an  ensemble  of  252  trajectories.  Time  zero  is  when  the  X3  is  first  released 
at  the  saddle  point  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  xenon  solvent 
while  the  dots  indicate  the  average  energy  for  a  reaction  in  the  absence  of  solvent  (gas  phase).  The 


translational  energy  arises  from  and  decays  to  the  solvent  in  0.15  ps,  and,  although  the  figure  does  not 
show  this,  vibrational  energy  starts  to  decay  in  >  100  ps. 


Figure  14.  Plot  of  the  average  energy  versus  time  for  the  reaction  X  +  X2  -»  X2  +  X  with  the  5 
kcal/mol  energy  barrier,  for  an  ensemble  of  64  trajectories.  Time  zero  is  when  the  X3  is  first  released  at 
the  saddle  point  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  helium  solvent 
while  die  (tots  indicate  the  average  energy  for  a  reaction  in  the  absence  of  solvent  (gas  phase).  The 
translational  and  rotational  energy  arises  from  and  decays  to  the  solvent  over  a  time  scale  measuring 
longer  than  1.0  ps,  and,  although  the  figure  does  not  show  this,  vibrational  starts  to  transform  in  less 
than  10  ps. 


Figure  15.  Plot  of  the  average  energy  of  versus  time  for  the  reaction  X  +  X2  -*  X2  +  X  with  the  5 
kcal/mol  energy  barrier,  for  an  ensemble  of  252  trajectories.  Time  zero  is  when  the  X3  is  first  released 
at  the  saddle  point  The  solid  lines  show  the  average  energy  for  the  reaction  in  liquid  xenon  solvent 
while  the  dots  indicate  the  average  energy  for  a  reaction  in  the  absence  of  solvent  (gas  phase).  The 
translational  energy  arises  from  and  decays  to  the  solvent  in  0.15  ps,  and,  although  the  figure  does  not 
show  this,  vibrational  energy  starts  to  transform  in  >  100  ps. 


Figure  16.  Solvent  momentum  sphere  of  influence  on  reaction  with  the  20  kcal/mol  barrier  in  helium 
solvent  For  details,  see  the  caption  of  Fig.  5. 
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